{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap aggregation and link-confidence quantification with TIGRAMITE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TIGRAMITE is a time series analysis python module. It allows to reconstruct causal graphical models from discrete or continuously-valued time series based on the PCMCI framework and create high-quality plots of the results.\n",
    "\n",
    "This Nature Review Earth and Environment paper provides an overview of causal inference for time series in general: https://github.com/jakobrunge/tigramite/blob/master/tutorials/Runge_Causal_Inference_for_Time_Series_NREE.pdf\n",
    "\n",
    "This tutorial explains **Bootstrap aggregation (bagging) for time series causal discovery**, which is implemented as the function `PCMCIbase.run_bootstrap_of`. The bootstrap aggregation is a general meta-algorithm that can be combined with most causal discovery methods of TIGRAMITE, e.g., `run_pcmci`, `run_pcalg_non_timeseries_data`, `run_pcmciplus`, `run_lpcmci`, among others, including the whole range of conditional independence tests. You can refer to the preprint below for additional information.\n",
    "\n",
    "---\n",
    "**Preprint on bootstrap aggregation and confidence measure for time series causal discovery:**\n",
    "Debeire, K., Runge, J., Gerhardus, A., Eyring, V. (2023). Bootstrap aggregation and confidence measures to improve time series causal discovery. \n",
    "[https://arxiv.org/pdf/2306.08946](https://arxiv.org/pdf/2306.08946)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline     \n",
    "\n",
    "import tigramite\n",
    "from tigramite import data_processing as pp\n",
    "from tigramite.toymodels import structural_causal_processes as toys\n",
    "from tigramite import plotting as tp\n",
    "from tigramite.lpcmci import LPCMCI\n",
    "from tigramite.pcmci import PCMCI\n",
    "from tigramite.independence_tests.parcorr import ParCorr\n",
    "from tigramite.independence_tests.cmiknn import CMIknn\n",
    "from tigramite.pcmci_base import PCMCIbase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrap aggregation (bagging) and confidence estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In bootstrap-based bagging, a random sample in the training set is selected with replacement\n",
    "—meaning that in each bootstrap sample each data point can be drawn more than once. Several data samples are generated in this fashion to produce a set of replicates (also called resamples). The machine learning models are\n",
    "then trained independently on each replicate and, finally, the outputs are averaged for prediction tasks\n",
    "or aggregated for classification tasks (for example by majority voting). In our case, the training set is the input time series, the machine learning model is the causal discovery method and the output is a causal graph.\n",
    "\n",
    "The main interest of the bagging approach in time series causal discovery is to improve the robustness of the output graph and to provide confidence estimates for links in the graph. Since time series causal discovery is sensitive to temporal dependencies, it is essential to retain the temporal dependencies in the resampling procedure. However, standard resampling inevitably destroys (at least parts of) the temporal lagged dependencies. To resolve this problem, we employ the resampling strategy as illustrated below:\n",
    "\n",
    "<img src=\"figures/bootstrap_for_timeseries.PNG\" width=600px/> \n",
    "\n",
    "In the end, our bagging approach combined with a time series causal discovery method (here PCMCI+) can be summarized as follows: The original timeseries is resampled $B$ (number of bootstrap replicas) times. On each resample, PCMCI+ is run independently, resulting in one output graph. The $B$ output graphs are then aggregated by relative majority voting of the link types for each pair individually. Ties are resolved by the preference order (no link, ×−× and ◦−◦); and in case of a tie between → and ←, a conflicting link ×−× is returned.\n",
    "\n",
    "The returned results of Bagged-PCMCI+ are: \n",
    "\n",
    "1. An ensemble of $B$ causal graphs from applying PCMCI+ to $B$ datasets that are obtained by resampling\n",
    "with replacement while preserving temporal dependencies,\n",
    "\n",
    "2. the aggregation of all these graphs by majority voting (at the level of each individual edge) to a final output graph, \n",
    "\n",
    "3. link frequencies for the final aggregated graph that provide a confidence measure for the links.\n",
    "\n",
    "Below we represent this measure of confidence via the width of the links in the output graph: the confidence measure (proportionnal to the arrow widths) for $X_2 \\rightarrow X_3$ at time lag 2 is larger than for $X_4 \\rightarrow X_1$ at time lag 1.\n",
    "\n",
    "<img src=\"figures/bagged_pcmciplus_summary.PNG\" width=600px/> \n",
    "\n",
    "The bagged version has one further parameter compared to the chosen causal discovery algorithm: the number of bootstrap replicas $B$. In practice, we recommend to use $B$ as large as possible. The implementation uses parallelization via ``joblib``. In our numerical experiments $B=25$ gave already good results, but we recommend using $B\\geq 100$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section demonstrates and explains the application of Bagged-PCMCI+ on synthetic data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data generation model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 640x480 with 1 Axes>, <Axes: >)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4AElEQVR4nO3deXhU5d3/8c9MFpKQQAj7EkCgxbXKpmVTkPUnmxAUWbQEKApFaVkEERRFZBF4BGrFIrKjgnFDAYsUkEUFEaECIsgihLAlIftkksz5/eEjj1SQzMyZzJzM+3VdvVqTub/3lyYyn7nPfe5jMwzDEAAACFp2fzcAAAD8izAAAECQIwwAABDkCAMAAAQ5wgAAAEGOMAAAQJAjDAAAEOQIAwAABDnCAAAAQY4wAABAkCMMAAAQ5AgDAAAEOcIAAABBjjAAAECQIwwAABDkCAMAAAQ5wgAAAEGOMAAAQJAjDAAAEOQIAwAABDnCAAAAQY4wAABAkCMMAAAQ5AgDAAAEuVB/NwCgZLlcLuXn58vhcKioqEgul0uSZLfbFRISooiICJUpU0Z2O58VgGBBGABKOcMwlJ2drYyMDOXl5Sk/P79Y48qUKaPIyEiVL19e0dHRstlsPu4UgL/YDMMw/N0EAPMVFBQoPT1daWlpKiws9KpWaGio4uLiVKFCBYWFhZnUIYBAQRgAShmXy6Vz584pNTXVJ/UrVqyoqlWrchkBKEUIA0Apkp2drdOnT3u9EnA9oaGhqlWrlqKjo306D4CSQRgASgGXy6UzZ87o0qVLJTpvbGysatSowSoBYHGEAcDiCgsLdfLkSeXl5fll/sjISNWtW1chISF+mR+A9wgDgIUVFhbq+PHjxb5DwFfKlCmjG264QaGh3KAEWBFre4BFFRUV6cSJE34PApKUn5+vEydOqKioyN+tAPAAYQCwqJSUFDkcDn+3cZnD4VBKSoq/2wDgAcIAYEGZmZklvlmwOC5duqTMzEx/twHATYQBwGIKCwuVnJzs7zauKTk52ee3NgIwF2EAsJiUlJSAvjZfVFTE5QLAYggDgIU4nU5lZGT4u43rysjIkNPp9HcbAIqJMABYSFpamr9bKLb09HR/twCgmAgDgEW4XC5LhYHU1NTLj0cGENgIA4BFZGZmlsib69mzZ5WYmKgePXqoV69e+uSTTzyq43K5uLMAsAhOIAQsIjk5uUSW3i9cuKDU1FTdeOONSk1N1YMPPqi1a9cqKirK7VpxcXGqUaOGD7oEYCbODgUsIjc3t0TmqVy5sipXrizpp8cVly9fXpmZmR6FgZLqGYB3uEwAWIDL5fL62GGXy6Vu3bppzpw5V3x9x44datSo0VUvBxw4cECGYahatWoezelwONg3AFgAYQCwADOOHbbb7RoyZIhWr159+fbEw4cPa/To0Ro5cqQ6dep0xesvXbqkCRMm6Nlnn/Vq3kA6MhnA1REGAAsw60S/Ll26KDY2VqtWrdLZs2c1fPhwde3aVQMHDrzidU6nUyNHjtSQIUN0xx13eDUnpxECgY89A4AFmLXUHhoaqkGDBmnevHnauHGjbr75Zj311FNXvMYwDE2cOFF33XWXunXr5vWc7FEGAh8rA0CQ6dq1qxwOhwzD0IwZMxQSEnLF9/fu3asNGzbo3//+t3r37q3evXvr+++/93g+wgAQ+FgZACzAZrOZVmvq1KmSftoT8N9BQJIaN26s/fv3mzaf3c5nDiDQ8W8pYAFXe9P2xPz587Vt2zatXLlSRUVFevfdd02p+1vM6h2A7xAGAAuIjIz0ukZSUpKWLVum+fPnq2HDhhowYIAWL16sgoICEzq8toiICJ/WB+A9wgBgASEhIQoLC/N4/LZt2zR16lRNmzZNt99+uySpX79+ys7O1kcffWRWm78SFhbGygBgAYQBwCI8OQFQ+ungoNGjR2vUqFFq37795a9HR0erX79+WrRokYqKisxq8wqe9gygZPFsAsAiUlNTlZKS4u823FK9enVVrFjR320AuA5WBgCLKF++vL9bcJsVewaCEWEAsIjQ0FDFxsb6u41ii42NVWgody8DVkAYACzESkvuVuoVCHaEAcBCIiMjTbnN0Nes0ieAnxAGAIupWbOmv1u4Liv0COD/EAYAi4mIiFDVqlX93cY1Va1alYOGAIshDAAWVKlSpYBcho+MjFSlSpX83QYANxEGAAuy2WyKj48PqNP9QkJCFB8fb+pDlQCUDMIAYFHh4eEqKipSRkaGv1uR3W7XDTfcoPDwcH+3AsADhAHAorKzs/XAAw/ogQce0O7du/3WR0hIiOrVq8c+AcDCCAOABRmGoWHDhum7775TSkqKBg8erPnz5/v8CYT/rUyZMqpfvz5BALA4nk0AWNAbb7yhwYMH/+rrjRo10tSpUxUfH+/zHqpUqaJKlSrJbuczBWB1hAHAYr799lvdeeedysvLu+r3Y2Ji9Oyzz6pTp04+mT8yMlI1a9ZkNQAoRQgDgIVkZ2erWbNm+u6776772j59+ujpp582bXd/dHS04uLiFBMTwx0DQClDGAAs5E9/+pOWLVtW7NfXq1dPy5YtU4MGDZSRkSGn0+nWfOHh4SpfvrwqVKjAnQJAKcYjxQCLWLJkiVtBQJKOHTumjRs3qmXLlqpatapcLpccDofy8vLkcDhUVFQkl8sl6afbA0NCQhQREaHIyEhFRESwHwAIEqwMABZw4MABNWvW7Jr7BK6lbdu22rhxY0AdTgQg8BAGgACXk5OjZs2a6dChQ26Nq1Klir755htVr17dR50BKC1YAwQCmGEYGjJkiNtBwGazaeXKlQQBAMVCGAAC2AsvvKC33nrL7XETJ05U+/btfdARgNKIywRAgEpKSlLv3r3dHnfPPfdo06ZN7BMAUGyEASAAff3112rVqpXbGwYrV66sb775RjVq1PBRZwBKIy4TAAEmJSVF3bt3dzsI2Gw2rVixgiAAwG2EASCA5OXl6f7771dycrLbYydMmKCOHTv6oCsApR2XCYAAYRiG+vfvrzfffNPtsXfffbc2bdqk0FDOEQPgPlYGgADx4osvehQE4uPjtXr1aoIAAI+xMgAEgHfffVcJCQluj4uKitKOHTt0xx13mN8UgKDBygDgZ3v37tXDDz/s0dgVK1YQBAB4jTAA+NHZs2fVvXt35ebmuj126tSp6tmzpw+6AhBsuEwA+InD4VCbNm305Zdfuj22f//+Wr58uWw2mw86AxBsWBkA/ODnZw54EgTuuusuvf766wQBAKYhDAB+MH36dK1cudLtcfHx8Xr//fcVERHhg64ABCsuEwAl7P333/foWj93DgDwFVYGgBK0b98+DRgwwKOxy5cvJwgA8AnCAFBCzp49q27duiknJ8ftsS+88IJ69erlg64AgMsEQIlwOBxq27atvvjiC7fH9uvXTytWrGDDIACfIQwAPlZUVKR+/fpp9erVbo+98847tWXLFkVGRvqgMwD4CWEA8CGXy6VBgwZp6dKlbo+tVauWdu3aperVq/ugMwD4P+wZAHzEMAyNGDHCoyAQFRWlDz/8kCAAoEQQBgAfMAxDY8aM0auvvurR+OXLl6tRo0YmdwUAV0cYAHzg2Wef1Zw5czwaO2XKFO4cAFCi2DMAmGzatGmaMGGCR2P79u2rlStXcucAgBJFGABMNHfuXP31r3/1aGyzZs20detW7hwAUOIIA4BJFi5cqKFDh3o0ljsHAPgTewYAE6xYsUKPPvqoR2OrVq2qTZs2EQQA+A0rA4CX3nnnHfXp00cul8vtsRUrVtSWLVt06623+qAzACgewgDghY8//lj333+/CgsL3R5bvnx5/fvf/1bjxo190BkAFB+XCQAPffrpp0pISPAoCJQtW1br168nCAAICIQBwAPbt29Xjx49lJ+f7/bYiIgIffTRR2revLkPOgMA9xEGADft3r1b9913n3Jzc90eGx4ervfff19t2rQxvzEA8BBhAHDDvn371KlTJ2VlZbk9NiQkRKtXr1anTp180BkAeI4wABTToUOH1KFDB6Wnp7s91m63a+XKlerRo4cPOgMA7xAGgGL44Ycf1L59e124cMGj8YsWLVKfPn1M7goAzEEYAK7jP//5j+6++26dOXPGo/GvvPKKBg4caG5TAGAiwgDwG7Zu3arWrVt7HARmzZql4cOHm9wVAJiLMABcQ1JSkjp16qSMjAyPxj///PMaPXq0yV0BgPkIA8BVvPrqq3rggQc8OkdAksaNG6eJEyea3BUA+AZhAPgFwzA0adIkDR8+XJ6e1P34449r2rRpstlsJncHAL7BswmA/1VYWKjHHntMixYt8rjGkCFD9Nprr8luJ2cDsA7CACApNzdXDz30kNauXetxjf79+2vp0qUKCQkxsTMA8D3CAIJeamqqunXrps8//9zjGgkJCXrrrbcUGhpqYmcAUDJYy0RQ+/HHH9W6dWuvgsCAAQO0atUqggAAyyIMIGj95z//UfPmzXXo0CGPa4wdO1ZLly5VeHi4iZ0BQMniowyC0meffabu3bt7fIaAJM2ZM0d/+9vfTOwKAPyDMICg895776lv374enyEQFhampUuXqm/fviZ3BgD+wWUCBJUFCxaod+/eHgeB6OhorVu3jiAAoFQhDCAoGIahZ555RsOGDZPL5fKoRpUqVbR161a1b9/e5O4AwL+4TIBSr7CwUMOHD9fChQs9rtGgQQNt2LBB9evXN7EzAAgMhAGUapcuXdKAAQP08ccfe1yjSZMmWrdunapUqWJiZwAQOLhMgFLrm2++UZMmTbwKAh07dtSWLVsIAgBKNcIASqUlS5aoefPmOnbsmMc1+vfvr7Vr1yo6OtrEzgAg8BAGUKo4HA4NHTpUiYmJcjgcHtcZM2aMli1bxmFCAIICewZQapw4cUK9e/fWnj17vKoze/ZsjRo1yqSuACDwEQZQKqxfv179+/dXenq6xzXCwsK0ZMkS9evXz8TOACDwcZkAluZyuTR58mR16dLFqyAQHR2tjz/+mCAAICixMgDLSk1N1YABA7Rhwwav6lSpUkXr1q1TkyZNTOoMAKyFMABL2r17t3r37q0ff/zRqzr169fXJ598wmFCAIIalwlgKYZh6J///KdatWrldRC47777tGvXLoIAgKDHygAsIzc3V8OHD9fSpUuv+Hp0dLSaN2+uKlWq6NSpU9q5c6cKCwuvWcdms+m5557T008/LbudPAwAhAFYwtGjR5WQkKD9+/df8fXOnTvrvffeU2Zmpk6dOqXf//73SklJUc+ePXXw4MFf1alYsaJWrVqljh07llTrABDwbIZhGP5uAvgtH3zwgf70pz8pIyPjV99r1qyZoqOjtXnzZklShQoVtHnzZiUnJ6tLly6/eu2aNWtUp06dEukbAKyCNVIErMLCQj311FO6//77rxoEpJ82Ev4cBCQpPT1dq1evVrNmza543bBhw7Rt2zaCAABcBZcJEJCOHz+uxMREbd261e2xTZo0uXyJoEyZMlq4cKEefvhhs1sEgFKDMICA4nK59Oqrr2rcuHHKyclxe3zv3r3Vo0cP3XPPPQoLC9MXX3yhO+64w/xGAaAUIQwgYPzwww8aPHiwR6sBknTvvfdq2bJlGjZsmHbs2KH169cTBACgGNgzAL9zuVyaN2+e/vCHP3gcBO655x6tXbtW48eP18KFC7VkyRJ17tzZ5E4BoHRiZQB+deTIEQ0aNEjbt2/3uEarVq308ccfa/LkyZo3b54kaevWrcrNzdWjjz7KWQIAcB3cWgi/KCoq0rx58zRhwgQ5HA6P6zRq1Ehbt27Vpk2bNHPmzCu+9/nnn+vTTz9Vu3btvG0XAEo1VgZQ4g4fPqzExER9/vnnXtdq1qyZnE6nWrdurdatW1/xvUqVKunEiRNezwEApR0rAygxRUVFmjNnjiZNmqT8/Hyv64WGhv7mscN2u10HDx5Uw4YNvZ4LAEozVgZQIg4ePKhBgwbpyy+/NKXejTfeqDVr1ujixYu6dOnSVV/TuHFj1a5d25T5AKA0IwzApwoLC/XSSy9p8uTJcjqdptRMTEzU3LlzFRMTY0o9AAh2hAH4zLfffqvExER99dVXptSrVauWFi5cyC2DAGAy7rmC6QoKCvTCCy+ocePGpgWBIUOG6NtvvyUIAIAPsDIAU+3fv18DBw7U3r17TalXu3Ztvf766+rQoYMp9QAAv8bKAExx/vx5DR8+XI0bNzYtCDz22GP69ttvCQIA4GOsDMArDodDc+fO1YsvvqjMzExTatatW1eLFi3Svffea0o9AMBvIwzAI4ZhaPXq1Ro3bpxOnjxpWt0RI0Zo2rRpio6ONq0mAOC3EQbgts8//1yjRo3SF198YVrN+vXra9GiRbrnnntMqwkAKB72DKDYjh8/roceekgtWrQwLQjYbDb99a9/1b59+wgCAOAnrAzgujIyMvTiiy/q5ZdfNu3gIEn63e9+p8WLF6tly5am1QQAuI+VAVxTYWGh/vGPf6hBgwaaOXOmaUHAZrNp9OjR+uabbwgCABAAWBnArxiGoXXr1mns2LE6dOiQqbVvvPFGvfHGG2revLmpdQEAnmNlAFfYt2+fOnbsqK5du5oaBOx2u8aNG6e9e/cSBAAgwLAyAElSSkqKJk2apDfeeENmP9X6zjvv1Pz583XnnXeaWhcAYA7CQJDLzc3V7NmzNWPGDOXk5Jhau3bt2po+fbr69Okju51FKAAIVISBIJWdna2FCxdq1qxZOnPmjKm1Y2JiNGHCBI0cOVKRkZGm1gYAmI8wEGQuXryo+fPna/78+UpPTze1tt1u15///Gc999xzqlq1qqm1AQC+QxgIEidPntScOXO0cOFC5eXlmV6/c+fOeumll3TrrbeaXhsA4FuEgVLuwIEDmjlzplatWqXCwkLT699yyy2aPXu2OnXqZHptAEDJIAyUUp9//rmmT5+uDz/80Cf1q1SpoilTpmjQoEEKDeXXCACsjL/FSxHDMLRhwwZNnz5dn332mU/miIiI0KhRozR+/HjFxMT4ZA4AQMkiDJQChYWFWrNmjaZPn679+/f7bJ7+/fvrxRdfVO3atX02BwCg5BEGLCwvL09LlizRSy+9pOPHj/tsnlatWmn27NkcGgQApRRhwIIuXbqkV199VS+//LLOnz/vs3nq1aunmTNnqlevXrLZbD6bBwDgXxwL5wcOh0Nvv/222+NSUlI0btw41a5dWxMmTPBZEChfvrxmz56tgwcPKiEhgSAAAKWczTD7IHo/KiwslMPhUGFhoVwul6SfHpdrt9sVERGh8PBwv7+xFRUVqU+fPtq6dauSk5MVHh7+m683DEOfffaZXnvtNSUlJZn2GOGrqVy5skaOHKnhw4erQoUKPpsHAKzAMAw5nU45HA65XK7Lz22x2+0KDQ1VREREqbmbytJ/CqfTqYyMDOXm5iovL++699HbbDZFREQoMjJSMTExio6OLtFwYBiGnnjiCSUlJUmS1q5dq4SEhKu+Nj09XcuWLdOCBQv03Xff+bSvOnXqaOzYsUpMTFRUVJRP5wKAQGUYhrKzs5WVlaW8vDw5HI7rPrgtNDRUkZGRioqKUvny5a/7AS9QWW5lwDAMZWVlKS0tTdnZ2V7VCgsLU1xcnGJjYxUWFmZSh9c2depUTZw48fI//7//9/+0bt26y/9sGIa+/PJLLViwQG+//bYcDodP+7n11ls1fvx4PfjggyXy5weAQFRQUKD09HSlpaV5fThbdHS04uLiFBMT4/eVaHdYJgwYhqHMzEylpKT45CS9ChUqqFq1agoJCTG9tiQtWrRIQ4YMueJrdrtdJ0+eVPny5bVy5UotWLBA+/bt88n8v9SqVSuNHz9e9913n6V+WQHATEVFRTp79qzpz2mRfloxqF69usqVK2eJv2ctEQYKCgp05swZZWVl+XSekJAQ1axZU+XKlTO17ttvv61+/fpd3sfwS02bNtV3333n9SpHcXTt2lXjxo1Tq1atfD4XAASyjIwMnTlzRkVFRT6dJyYmRjVq1Aj41deADwPp6ek6c+bMda/bmKlcuXKqUaOGKRtDPvjgA/Xu3dsnqxnFERISor59++rJJ5/Ubbfd5pceACBQFBYW6syZM8rMzCyxOW02m2rUqBHQG7MDNgwYhqHz58/rwoULfpk/PDxcdevW9WozyCeffKLu3bv79A6Aa4mIiNCQIUM0evRo1a1bt8TnB4BA43Q6deLECb/8nSz9dMdWlSpVAvKyQUCGAcMwlJKSorS0NL/2ERoaqnr16nkUCLZu3arOnTv7fBPgf4uNjdWIESP0+OOPq0qVKiU6NwAEKqfTqWPHjvltlfZncXFxql69esAFgoAMAykpKUpNTfV3G5J+uuOgXr16bl3v+eKLL9ShQ4cS2Qfwsxo1amjUqFEaOnQoDxACgF8oKCjQsWPHVFBQ4O9WJEkVK1ZU9erV/d3GFQIuDFy6dEmnT5/2dxtXiIqK0g033FCsJLd37161bdtWGRkZJdCZ9Pvf/15PPvmkBgwYoDJlypTInABgFYZh6NixY8rLy/N3K1eIj49X+fLl/d3GZQF16NDPdw0EmtzcXKWmpqpSpUq/+boDBw6oY8eOJRIEOnXqpMcee0zdunXz2e2QAGB1qampARcEJCk5OVlRUVEBc5dBwIQBwzB0+vTpq95+FwjOnj2r6OhoRUREXPX7R44cUfv27XXx4kWf9VC5cmUNHjxYf/7zn1WvXj2fzQMApYHD4dDZs2f93cZVuVwunT59WnXr1g2I/QMBEwYuXbqknJwcf7fxm06dOqUGDRr86gd3/PhxtWvXzme/dG3atNFjjz2mnj17WvaoSwD4pZkzZ6pGjRp66KGHfHK+v2EYOnXqlOl1zZSTk6NLly4FxC2HAfHUQsMw/HYLoTvy8/N/FViOHj2qe+65x/RfugoVKuhvf/ubDh06pM2bN6tPnz4EAQClxtGjR/Xwww/rlltu0YoVK0zf5Z+Tk6P8/HxTa/rChQsXSvQcnWsJiDCQnZ3tt/s+3fXLywDff/+96UHgpptu0tKlS5WcnKw5c+boxhtvNK02AASa77//3iehwJeXbM3kdDoDYlU8IMJAoNxGWBw/B5dDhw7pnnvuMX3DY506dfTII48oMjLS1LoAEMjMDAX5+fklemu3twIhuPg9DDidzhL5oeXk5Oihhx5S79691bNnT73zzjse1/ryyy/Vpk0bn+wR+OSTTwL+OhcA+IoZocAXDx66lpEjR6pFixYaNWqUxzUCYXXc7+cMpKenKzk52efzFBUVyel0KjIyUnl5eerZs6feeustxcbGul2rTJkyOnLkiHbt2qVz585d/s/58+d17tw5rx+olJCQoJ49e3pVAwAC2euvv64tW7Zc93W/+93vNGnSJPXt27fYGw2PHj1aYqe/7tq1S7m5ufrwww81Z84cj+vUqlXLo/cjs/j9boKSuv8zJCTk8tK70+mUy+XyeNNGfn6+Vq9eraVLl5rZ4mVJSUlKSkrySW0AsJIjR47okUce0ZQpU4oVClwuV4keA3/nnXdq9+7dXtfJy8vzaxjw+2UCb8OAy+VSt27dfpXIduzYoUaNGumTTz65/LXMzEwlJCSoffv2SkxM9Op2Dl8/9hIA8H9+DgU333yzli9ffs3LB2bcQeDO+4pZcnNzTa/pDr+GAcMwvE5wdrtdQ4YM0erVqy+f/Hf48GGNHj1aI0eOVKdOnS6/tly5ckpKStL69eu1bt06rzZtEAYAoORdLxSYsdrszvuKWRwOh19vMfRrGCgsLDTlD9+lSxfFxsZq1apVOnv2rIYPH66uXbtq4MCBV319pUqV9Pvf/1579uzxeM5APSkRAILBtUKBWQ8jcvd9xVuGYfj1Q6bfVwbMEBoaqkGDBmnlypUaPny4br75Zj311FNXvObixYuX71rIzs7Wnj17VLduXY/nDIRDIgAg2P13KDBrV35x3lfM5s8PmaUiDEhS165dLy+zzJgx41cP7zl37pwGDhyohIQEPfLII+rbt68aNmxo2vwAAP85cuSI/vnPf+rw4cOm1bze+4okPfrooxo9erS2bdumdu3a6dtvv/V4Pn9+yPTr3QR2u3lZZOrUqZJ+esbB1X5gt9xyi1dnC/y3QHiwBABAatWqlZ577jm1bdtW586dM+0Qn+u9r0jSa6+9ZspckrnviW7P7beZZd4ffP78+dq2bZtWrlypoqIivfvuu6bU/S2EAQDwr1atWmnTpk367LPPdO+998pms5n2SHd/vK8EbRgICQnx+lnOSUlJWrZsmebPn6+GDRtqwIABWrx4sWmbSK7FrF84AIB7rhYCfnatx8y7wx/vK2FhYX59X/H7OQPenMG/bds2TZ06VdOmTdPtt98uSerXr5+ys7P10UcfmdXiVREGAKBktWzZ8poh4GfePtfFX+8rUVFRPqtdHH4/jvjixYsenfF/4MABJSYm6oknntCAAQOu+N7f//53bdiwQR988IFP3rRDQ0P1/vvva/369abXBoDS7tChQx49tr58+fKqX7++GjRocPm/f/7f1atXv7zM/t1333n0TAN/vq9Uq1ZNlSpVMr1ucfk9DOTk5Oj48eP+bMFt5cqVU+3atf3dBgBY0qZNm9S+fXtTa0ZGRqpevXpq0KCBnnjiCVWpUsXU+r52ww03qGzZsn6b3++XCaKiorzeN1DSvDnGGACCXdu2bRUTE2Nqzby8PB04cEDNmzdX06ZNTa3ta2FhYX6/TOD3MGCz2VSxYkV/t1FsYWFhio6O9ncbAGBZDofDJ3/vz5kzR+PGjVNMTIylPmRWrFjR73eo+T0MSD990vb3/xHFFQg/NACwooMHD2rkyJGqWbOmTpw4YWrtefPm6W9/+5ska33ItNlsAbHa7PdHGEs/7cyvUKGC0tLS/N3KbwqUHxoAWEV+fr7effddLViwQJ999plP5vjHP/6hYcOGXfG1ChUq6Ny5cwF/dHyFChUC4u60gAgDklSlShVlZGQE9NMAq1WrFhA/NAAIdEePHtU///lPLV682LQTAa9m/vz5vwoC0k8fMqtVq6aUlBSfze2tkJCQgNnoGDBhIDQ0VDVr1tSPP/7o71auqmzZsoqLi/N3GwAQsAoKCrR27VotWLBAGzdu9Pl8M2bM0IgRI675/bi4OGVmZionJ8fnvXiiZs2aCg0NjLfhwOjif5UrV06xsbG6dOmSv1u5gt1uV61atdgrAABXkZmZqQULFmju3Lk6c+ZMicz5zDPP6Mknn/zN19hsNtWqVUtHjhwJuMfOx8bGqly5cv5u47KACgOSVL16deXk5Pj8OGF31KxZ01I7UwGgJJw7d05z587VP/7xD2VkZJTYvGPGjNHkyZOL9dqwsDDVrFlTp06d8m1TbggLC1P16tX93cYV/H7o0NUUFBTo2LFjAREIatSoweUBAPiFY8eOadasWXrjjTeUn59fonMPHz5cf//7391eqU1LSyuxVYvfEhYWpnr16gXcB8yADAOS5HQ6deLECTmdTr/1QBAAgP+zb98+zZgxQ2+//bZflt0HDhyoRYsWefx0P38HgvDwcNWtW1fh4eF+6+FaAjYMSFJhYaFOnjypvLy8Ep3XZrMpPj4+oK7nAIA/GIah7du3a/r06Vq3bp3f+ujTp49Wrlzp9R1dmZmZOnXqVInfchgZGak6deoEzIbB/xbQYUD66RcxNTW1xO4XLVu2rGrWrBmQyQ0ASorL5dLHH3+s6dOna+fOnX7tpUePHlqzZo1pS+tOp1PJycklcpeBzWZT1apVA/7AuoAPAz/z9Q/PbrerRo0aKl++fED/wADAlwoKCvTWW29pxowZOnDggM/na9SokYYNG6Zq1aqpe/fuv/p+hw4dtHbtWpUpU8bUeQ3DUEZGhs6cOeOzSx5W+nBpmTAg/fTDy8rKUlpamrKzs02pGRYWpri4OFWoUCFgl28AwNdyc3O1aNEizZo1y+fnvURGRqpfv3569NFH1bRpU9lsNrlcLv3ud7/TsWPHLr+uadOm2rx5s0+fB1NYWKj09HSlpaWZtmk9OjpacXFxiomJscyHS0uFgV9yOp1KT09Xenq628+tttlsiomJUVxcnMqWLWuZHxYAmC0tLU2vvPKK5s2b59OTAiXplltu0WOPPaYBAwYoNjb2V9+fOnWqJk6cKElq0KCBduzYUWIn9BmGoZycHKWlpSkrK8vty9KhoaGqUKGCKlSoYImVgP9m2TDwS4WFhXI4HMrLy1NeXp4KCwvlcrlkGIbsdrvsdrsiIiIUGRmpiIgIlSlThgAAIKidPn1a//M//6PXXnvNp9fOy5QpowceeECPPvqoWrZs+Zt/954+fVp16tRR5cqVtXPnTtWrV89nff0WwzCUn59/+X3F4XDI5XLJ5XLJZrPJbrcrNDRUkZGRl99XrL6yXCrCAACgeFJSUvTss89qyZIlPj3LpVq1avrrX/+qwYMHq1KlSsUe99BDD2ncuHFq1KiRz3rDrxEGACAI5Obmavbs2ZoxY4ZPVwLq16+vJ598Uo888ogiIiLcHl9QUBBwB/IEA8IAAJRiLpdLK1as0IQJE5ScnOyzeRo1aqTx48crISGBp7takLUvcgAArmnLli0aPXq0vv76a5/N0bZtW40fP14dOnRgL5aFEQYAoJT5/vvv9eSTT+qDDz7w2Rw9e/bUuHHjdNddd/lsDpQcwgAAlBJpaWl6/vnn9corr7h9y3VxhIaG6uGHH9bYsWN10003mV4f/kMYAACLczqdeuWVV/T888/r0qVLptePiorS0KFDNWrUKMXHx5teH/5HGAAAizIMQ++9956efPJJ/fDDD6bXj4uL0xNPPKERI0aoYsWKptdH4CAMAIAF7d69W6NHj9a2bdtMrx0fH6/Ro0dryJAhKlu2rOn1EXgIAwBgIadOndKECRO0YsUK02vfdNNNGjdunPr27WvJI3XhOcIAAFhAVlaWZsyYodmzZ8vhcJhau169epo+fboSEhJkt9tNrQ1rIAwAQAArKirSG2+8oUmTJuncuXOm1o6NjdWkSZP0l7/8xfRHBMNaCAMAEKAOHDigxMRE7d6929S6oaGhGj58uJ555hk2BkKSxHoQAASYgoICTZ06VY0bNzY9CPTo0UMHDhzQ3LlzCQK4jJUBAAgg+/fvV2JioulHCDdq1Ehz5sxRmzZtTK2L0oGVAQAIAAUFBXr++efVtGlTU4NAjRo1tGTJEn311VcEAVwTKwMA4GfffPONBg4cqH379plWMyoqSuPGjdPo0aM5KwDXRRgAAD9xOp164YUXNG3aNNOeJWCz2ZSYmKgpU6aoRo0aptRE6UcYAAA/2LNnjxITE/Wf//zHtJrt2rXTrFmzdMcdd5hWE8GBPQMAUILy8/P19NNP66677jItCNx444366KOPtHHjRoIAPMLKAACUkF27dikxMVEHDx40pV6lSpU0efJkDR06VGFhYabURHBiZQAAfMzhcGjcuHFq3ry5KUHAZrNp5MiROnLkiP7yl78QBOA1VgYAwIc+//xzDRo0SN99950p9Ro2bKg33nhDLVq0MKUeILEyAAAeyc/PV25u7jW/n5eXpzFjxqhly5amBAG73a6xY8dq7969BAGYjjAAAMWUlJSkBx98UPXr11dERITKli2rFi1ayOl0XvG67du36/bbb9fs2bNlGIbX8950003auXOnZs6cqcjISK/rAf/NZpjxmwoApdzBgwd1yy23qEWLFvrjH/+o2267TTk5ORoxYoQ2bdqke++9Vy6XSy+88IImT55sSgiw2+0aN26cnnnmGUVERJjwpwCujj0DAFAMJ06ckCStWrVKqampCgkJUcWKFTVixAg5nU6lpqbq4Ycf1vr1602Z79Zbb9XixYvVtGlTU+oBv4XLBABQDDfccIMkqUGDBmrSpIlWrVp1+XtHjhxRkyZNTAkCISEhmjhxor766iuCAEoMKwMAUAw33XSTNm/erMOHD2vy5MlXfG/UqFGmHCf8hz/8QYsXL1bjxo29rgW4g5UBACimNm3a6NFHH1V0dPQVX/c2CISGhmry5MnavXs3QQB+wcoAALipoKDAtFp33HGHlixZottvv920moC7WBkAADesXbtWp06d8rpOWFiYpkyZol27dhEE4HeEAQAohqKiIj399NPq3r27XC6XV7WaNGmiPXv2aOLEiRwljIDAZQIAuI4LFy6ob9++2rRpk1d1wsPDNXnyZI0dO1ahofz1i8DBbyMA/IYvvvhCDzzwgE6fPi1Jeuqpp9SkSRM1bNhQFStW1DvvvCNJGj58uM6fP3/NOvXr11dSUhKXBBCQCAMAcBWGYeiVV17RqFGjrtgwePToUUnSnj17rnh9fn7+NWt1795dS5cuVWxsrE96BbzFccQA8F9ycnI0dOjQKw4W8oTdbteLL76osWPHym5nixYCFysDAPALhw8fVkJCgg4cOOBVncqVK+utt97Svffea1JngO8QVQHgfyUlJalZs2ZeB4HmzZtr7969BAFYBmEAQNBzuVwaP368evfuraysLK9qPfHEE9qyZYtq1qxpUneA73GZAEBQczqdGjhwoN58802v6pQtW1avv/66HnroIZM6A0oOYQBA0MrKylKvXr306aefelWnYcOGSkpK0i233GJSZ0DJ4jIBgKB07tw5tWnTxusg0Lt3b+3evZsgAEsjDAAIOkePHlWLFi309ddfe1wjJCREc+bM0erVqxUTE2Nid0DJ4zIBgKDy1Vdf6b777tOFCxc8rlGtWjWtXr1arVu3NrEzwH9YGQAQNP71r3+pTZs2XgWBu+++W3v37iUIoFQhDAAICitXrlSXLl2Uk5PjcY0xY8Zo06ZNqlatmomdAf7HZQIApd7s2bM1ZswYj8fHxMRoyZIl6tWrl4ldAYGDMACg1HK5XBo7dqzmzJnjcY1atWppw4YN3C2AUo0wAKBUcjqdSkxM9OphQzfffLM2bNig+Ph4EzsDAg9hAECpk5WVpYSEBG3cuNHjGi1bttSHH36ouLg4EzsDAhMbCAGUKufOnVPbtm29CgLdu3fXxo0bCQIIGoQBAKXGDz/8oJYtW2rPnj0e1xg6dKiSkpIUGRlpYmdAYCMMACgV9uzZoxYtWuiHH37wuMazzz6rBQsWKDSUK6gILvzGA7C8jRs3qlevXsrOzvZovN1u1yuvvKLHHnvM5M4AayAMALC0pKQk9e3bVwUFBR6NL1OmjN5880317NnT5M4A67AZhmH4uwkA8MTHH3+s+++/X4WFhR6Nj42N1YcffsjRwgh6hAEAlvTpp5+qa9euys/P92h8zZo1tWHDBt16660mdwZYD2EAgOVs27ZNnTt3Vm5urkfjb7rpJm3YsEG1a9c2uTPAmribAICl7Nq1S126dPE4CDRv3lzbt28nCAC/QBgAYBn79u1T586dlZWV5dH4bt266dNPP+UwIeC/EAYAWMKhQ4fUoUMHpaenezR+8ODBevfddxUVFWVyZ4D1EQYABLyjR4+qXbt2unDhgkfjJ06cqIULF3KYEHAN/JsBIKD9+OOPateunVJSUjwaP3fuXD3xxBMmdwWULqwMAAhYKSkpateunX788UePxr/00ksEAaAYCAMAAtKFCxfUvn17HT161KPxzz33nMaMGWNyV0DpRBgAEHDS09PVsWNHHTx40KPx48aN06RJk0zuCii9OHQIQEDJzMxUhw4dtGvXLo/GP/7445o7d65sNpvJnQGlF2EAQMDIzc1V586dtW3bNo/GDxkyRK+99prsdhY9AXcQBgAEBIfDoe7du2vjxo0eje/Xr5+WLVumkJAQkzsDSj/CAAC/KygoUEJCgtauXevR+F69euntt9/mHAHAQ6ylAfCrwsJC9e/f3+MgcN999+nNN98kCABeIAwA8BvDMDRkyBCtWbPGo/H33nuv3nnnHYWHh5vcGRBcCAMA/GbKlClaunSpR2NbtmypDz74QJGRkSZ3BQQf9gwA8Is1a9bowQcf9Ghs06ZN9emnn6p8+fImdwUEJ8IAgBK3Z88etW7dWnl5eW6Pve2227R582ZVrFjRB50BwYkwAKBEpaSkqFmzZkpOTnZ7bMOGDbV161ZVrVrVB50BwYs9AwBKTF5ennr06OFREKhXr542bdpEEAB8gDAAoEQYhqFBgwZp9+7dbo+tVauWNm3apJo1a/qgMwCEAQAl4oUXXtBbb73l9riqVatq06ZNqlu3rvlNAZDEngEAJSApKUm9e/d2e1zZsmW1Y8cO3X777T7oCsDPCAMAfOrrr79Wq1at3L5zwGaz6b333lOPHj181BmAn3GZAIDPpKSkqHv37h7dQvjiiy8SBIASwsoAAJ/Iy8tTmzZttGvXLrfHPvzww1q6dKlsNpsPOgPw3wgDAExnGIb69++vN9980+2xf/zjH7V582ZFRET4oDMAV8NlAgCmmzp1qkdBID4+Xu+//z5BAChhrAwAMJWndw5ERUVpx44duuOOO8xvCsBvYmUAgGn27t2rRx55xKOxK1asIAgAfkIYAGCKn+8cyM3NdXvs1KlT1bNnTx90BaA4uEwAwGsOh0Nt2rTRl19+6fbY/v37a/ny5dw5APgRKwMAvGIYhgYPHuxRELjrrrv0+uuvEwQAPyMMAPDKzJkztWrVKrfH1apVizsHgADBZQIAHtu2bZvatGkjl8vl1rioqCht375djRo18lFnANzBygAAj1y4cEF9+/Z1OwhI0vLlywkCQAAhDABwm8vl0iOPPKLk5GS3x06ZMkW9evXyQVcAPEUYAOC2l156SRs2bHB7XN++ffX000/7oCMA3mDPAAC3bN++XW3atFFRUZFb4+68805t2bJFkZGRPuoMgKcIAwCK7eLFi2rUqJFOnz7t1rhatWpp165dql69uo86A+ANLhMAKBaXy6U//elPbgeB8PBwvffeewQBIIARBgAUy6xZs7Ru3Tq3x82ePVtNmzb1QUcAzMJlAiAIuVwuORwOFRUV6ee/Amw2m0JCQhQRESG7/crPCTt37tTdd9/t9j6BhIQErVmzhhMGgQAX6u8GAPie0+lUZmam8vLylJeXJ6fT+ZuvDw8PV2RkpCIjI1WuXDlt3LjR7SBQr149LVq0iCAAWAArA0ApZRiGsrKylJaWpuzsbK9qRUdHa+fOnXr88ceVlpZ23deHhYVp586dXB4ALII9A0ApYxiG0tPTdfjwYf34449eBwFJys7O1h/+8AclJSWpZcuW1309+wQAa2FlAChFnE6nkpOTlZOT47M5ioqKtHz5cr388stXvXTQq1cvvfPOO1weACyEMACUAoZhKDU1VefOnVNJ/Su9f/9+TZgwQSdPnrz8tRtuuEFff/21YmNjS6QHAOYgDAAWZxiGkpOTdenSpRKfOysrSzNmzNAHH3ygsLAw7dixQ82aNSvxPgB4hzAAWJjL5dKpU6eUlZXltx4Mw9DGjRsVGRmpYcOG+a0PAJ4jDAAWZRiGTp06pczMTH+3IkkqV66c4uPj2SsAWBB3EwAWdeHChYAJApKUmZmpCxcu+LsNAB4gDAAWlJeXp/Pnz/u7jV85f/688vLy/N0GADcRBgCL+XmfQKA6deqUXC6Xv9sA4AbCAGAx58+fv+5xwv7kdDoDctUCwLURBgALKSoqUmpqqr/buK7U1FS3n2UAwH8IA4CFpKenl9ihQt74+UhkANZAGAAswjAMXbx40d9tFFtqaqolggsAwgBgGVlZWSosLCyx+fLy8tSxY0fNmjXLo/EFBQWmPCQJgO8RBgCL8OXDh65m4cKFuu2227yqQRgArIEwAFhEbm5uic118uRJHT9+XK1bt/aqTkn2DMBzhAHAAgzDkMPh8KqGy+VSt27dNGfOnCu+vmPHDjVq1EiffPLJ5a/NmjVLI0eO9Go+SXI4HOwbACyAMABYQH5+vtdvqna7XUOGDNHq1auVkZEhSTp8+LBGjx6tkSNHqlOnTpKkf//736pbt67q1q3rbdsyDEP5+fle1wHgW6H+bgDA9Zl1yFCXLl306quvatWqVerZs6eGDx+url27auDAgZdfs3//fq1fv17/+te/lJubq8LCQpUtW9bjJxI6nU5FRESY0j8A3+CphYAFZGRkmHYE8erVqzVv3jxVqVJFNWvW1Msvv6yQkJCrvvb999/X0aNHNWbMGI/ni4+PV/ny5T0eD8D3uEwAWICZmb1r166Xr+XPmDHjmkHALHzeAAIflwkAC7DZbKbVmjp1qiTp0qVL1w0C999/v9fzmdk7AN9gZQCwALM+vc+fP1/btm3TypUrVVRUpHfffdeUur/F1ysPALxHGAAswIwNeElJSVq2bJnmz5+vhg0basCAAVq8eLEKCgpM6PDa2DwIBD7CAGABoaGhCg31/Kretm3bNHXqVE2bNk233367JKlfv37Kzs7WRx99ZFabv+Jt3wBKBmEAsIjIyEiPxh04cECjR4/WqFGj1L59+8tfj46OVr9+/bRo0SKfPW44KirKJ3UBmItbCwGLuHDhgs6dO+fvNtxStWpVVa5c2d9tALgOVgYAi7DivfpW7BkIRoQBwCLCw8MVHR3t7zaKLTo6WuHh4f5uA0AxEAYAC6lUqZK/Wyg2K/UKBDvCAGAhZcuWtcSn7fDwcJUtW9bfbQAoJsIAYCE2m03Vq1f3dxvXVaNGDU4eBCyEMABYTExMjOLi4vzdxjXFxcVZam8DAMIAYEnVqlVTWFiYv9v4lbCwMFWrVs3fbQBwE2EAsCC73a74+Hh/t/Er8fHxstv5awWwGv6tBSwqKipKtWvX9ncbl9WpU4cTBwGLIgwAFlauXDnVrl3br5v1bDabateurZiYGL/1AMA7HEcMlAI5OTk6efKkXC5Xic5rt9tVp04dbiMELI4wAJQSBQUFSklJUWZmZonMV65cOdWoUYOnEgKlAGEAKGUyMzOVnJzssycRhoSEqGbNmipXrpxP6gMoeYQBoBQqKipSWlqaUlNTVVhYaErNsLAwxcXFKS4uTiEhIabUBBAYCANAKWYYhrKzs5WWlqasrCyPavx8yFF0dDSnCgKlFGEACBIul0sOh0MOh0N5eXnKy8tTUVGRfv4rwGazKSQkRJGRkYqMjFRERIQiIiI4NwAIAoQBAACCHJEfAIAgRxgAACDIEQYAAAhyhAEAAIIcYQAAgCBHGAAAIMgRBgAACHKEAQAAghxhAACAIEcYAAAgyBEGAAAIcoQBAACCHGEAAIAgRxgAACDIEQYAAAhyhAEAAIIcYQAAgCBHGAAAIMgRBgAACHKEAQAAghxhAACAIEcYAAAgyBEGAAAIcoQBAACCHGEAAIAg9/8BV+iihSEI7s0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set seed for reproducibility\n",
    "seed = 1111\n",
    "# Choose the time series length\n",
    "T = 100\n",
    "\n",
    "# Specify the model (note that here, unlike in the typed equations, variables\n",
    "# are indexed starting from 0)\n",
    "def lin(x): return x\n",
    "\n",
    "links = {0: [((0, -1), 0.3, lin), ((2, 0), 0.5, lin), ((3, -1), -0.5, lin)],   # X1\n",
    "         1: [((1, -1), 0.3, lin)],                                             # X2\n",
    "         2: [((2, -1), 0.3, lin), ((1, -2), 0.4, lin)],                        # X3\n",
    "         3: [((3, -1), 0.3, lin)]                                              # X4                            \n",
    "        }\n",
    "\n",
    "var_names = [r'$X^{%d}$' % j for j in range(1, len(links)+1)]\n",
    "\n",
    "# Show ground truth causal graph\n",
    "tp.plot_graph(\n",
    "    graph = PCMCI.get_graph_from_dict(links),\n",
    "    var_names=var_names,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now generate data from a structural causal process with this graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data according to the full structural causal process\n",
    "data, nonstationarity_indicator = toys.structural_causal_process(\n",
    "    links=links, T=T, seed=seed)\n",
    "assert not nonstationarity_indicator\n",
    "\n",
    "# Number of variables\n",
    "N = data.shape[1]\n",
    "\n",
    "# Initialize dataframe object, specify variable names\n",
    "dataframe = pp.DataFrame(data, var_names=var_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard PCMCI+\n",
    "\n",
    "We start by illustrating results with the standard PCMCI$^+$ method. Please consult the [tigramite overview tutorial](https://github.com/jakobrunge/tigramite/blob/master/tutorials/causal_discovery/tigramite_tutorial_causal_discovery_overview.ipynb) for selecting hyperparameters in practice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAHBCAYAAAD0E7h1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzBklEQVR4nO3dd5gc1Z3v/09Vp9FETVBGAeUwkkACBEKAMBbCIKLB2EJrwhovNuZZDCz32We518A167X9w4v9W+N7FxvBrg14CV7jJIKwBBLZBCExGEkgBMqT83R31bl/jGbQoDTT06Gq6/16Hj1qdXWf/mp6quvTp845ZRljjAAAQGDZuS4AAADkFmEAAICAIwwAABBwhAEAAAKOMAAAQMARBgAACDjCAAAAAUcYAAAg4AgDAAAEHGEAAICAIwwAABBwhAEAAAKOMAAAQMARBgAACDjCAAAAAUcYAAAg4AgDAAAEHGEAAICAIwwAABBwhAEAAAKOMAAAQMARBgAACDjCAAAAAUcYAAAg4MK5LgBAdhljuv+4rkz3Hd0bLEuWJMu2ZVmWLMvKYZUAsokwAOQ5Y4xc15XrOHJdV6bn4H8UlmXJtm3ZoZDs/QEBQH6yTH8/GQD4ijFGTjIpx3H6HQAOx7IshUIhhcJhQgGQhwgDQJ4xxiiZTMpJJjPSfigcVphQAOQVwgCQRxzHUSIez8prRaJRhUKhrLwWgMwiDAB5wBijRCIh13Gy+rp2KKRIJEIvAeBzhAHA54wxind1DXpcQKosy1I0FiMQAD7GOgOAj+U6CHilBgCDQxgAfMpLB2Ev1QJg4AgDgE8lEglPHXx7xi0A8B/CAOBDjuNkfbBgf7iOI8eDdQE4MsIA4DPGmKxNH0xFIh73VI8FgKMjDAA+4+Ug0MMPNQL4FGEA8JGe6wx43UCugQAg9wgDgI8kM7TEcCb4qVYg6AgDgE/0XHjIL5xkkt4BwCcIA4BPZGv2wMeffKKzly7V8fPm6cSTTtLjTzyRcltenPEA4GAsRwz4RCIez8q0vV27dmnv3r2aO3eu9u7dq1MWLtSGt99WUVHRgNsKhUKKRKMZqBJAOoVzXQCA/snWwMFRo0Zp1KhRkqThw4eroqJC9Q0NKYUBPwx2BMBpAsAXjDGDPv/uuq7mHnec/um22/rc/8wzz6i0rOyQpwP+8pe/yHVdjT3mmJReMx11A8g8wgDgA+k4oNq2rX+45Rbdd999amhokCRt2LBBV6xYoTvuuENfvOSSPo+vq6vT1669Vv/2b/82qNclDADex5gBwAccx0nLQj7JZFKz58zRiiuu0JVXXqkzFi/W+cuW6Z577unzuK6uLp23bJmuufpqLV++fFCvGYlGFQqFBtUGgMwiDAA+4CSTabsI0M9//nN95/bbNXr0aI0fN06//vWv+xysjTG68qqrNHXKFN32mVMKqYhEIgqFGZ4EeBlhAPCBdIaB1tZWjR03TpMmTdLaNWsOGhi4/sUXtWTJEs2uru697xe/+IWqD/j3QBAGAO9jDwX8wLLS1tS3b7pJklRXW3vI7vtTFy5Ue1tb2l4vnbUDyAwGEAIBcscdd2jVqlVau2aNko6jBx54INclAfAAwgDgA7Y9+F115cqV+vFPfqLHHntMc+bM0beuv14/+td/Tdvph8NJR+0AMou9FPABy7IG1d3+1FNP6cZvf1v3/+IXWnDSSZKkb3zjG2ppadFDDz2UrjIPZlndtQPwNMIA4BOpfsN+4403dMWKFbrrrrt00UUX9d5fWlqqb1x3nf6/u+/O2DLH9AoA/sBsAsAnksmkkhnu0k+3cCSiMDMJAM8jtgM+4ceFe/xYMxBEhAHAJyzL8lW3u23bjBcAfMI/nywAFI5Ecl1Cv/mpViDoCAOAj/jl27bfejGAoGNvBXwmEo3muoSj8kONAD5FGAB8xrZtT4/QD4fD9AoAPsMeC/hQKBz25OkCy7K4KBHgQ4QBwIcsy/JkV3wkGvVkSAFwZIQBwKds21Y0Fst1Gb2isRinBwCfYs8FfMwrgYAgAPgbyxEDecB1XSXicWV7d+45XUEQAPyNMADkCWOMksmknGQyK68XCocV9uhARgADQxgA8kymewnoDQDyD2EAyEPGGLmuq2QyKeO6aWnT2r++gV9WQQTQf4QBIM+5rivHceQ4jjTQ3d2yFAqFFAqF6AkA8hhhAAgQY4yM68rd/7cxRj0fAJa6TwFYti17/9/0AADBQBgAACDg6PcDACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAL57qAdDLGyLiujDHd/5ZkSZJlyTrgDwAAR2OM6f0jYz49pkjdxxPbzptjiq/DgDFGTjIp13Xlum6/nmPZtmzLkh0Kyc6jNxIAMDjGmO7jiePIPeCL5dHYti3bthUKh317TLFMf/+3HtHzZvWEgMGwLEuhUMjXbyAAYHB6vlg6jtPvAHA4PaHAb182fRUGHMdRIh7PSNt2KKRIJOKrNw8AkDpjjBKJhFzHyUj7kWhUoVAoI22nmy/CgDFGiXh80D0B/eGnNw8AkJpMfrk8kG3bikSjnv+i6fkw4CSTSiQSWX1Nv7x5AICByeaXywNFIhGFwt4dpufZMGCMUTKZlJNM5uT1LctSNBYjEABAnjDGKN7VNehxAakKhcMKe3SMmifDgDFGyURCTobO4wxELBaTZbMcAwD4mXFddXV15boMhUIhhT04Ps2TR7nk/lGdXtCVwxQJABg8Y4wngoDUPVYhmaMe7yPxXBhwHCdnpwYOJx6PEwgAwIeMMYpnYaDgQDge+sLbw1NhoGdgh9cY1/XcGwcAODrHcWSyPFiwPxIe+5LpmTDg1SDQI5lIZH30KQAgda7rKpnl2WgD4aVA4Jkw0LP8o5d56Y0DABye179gSupd+tgLPBEGeqYRel3PUsgAAG8byLUFcimZTHqiTk+EAb+8aZI8N7gRAHAwv3xWe+VLpifCgF/eNEkDukIiACD7/PY57YVjYM7DgMnSm9bS0qJTFy3SggULdMKJJ+r+++9PuS1mFgCAd2XzM/pLl1+uUaNH6yvLl6fchuu6OZ/xkPMVCLN17QHHcdTV1aXCwkK1t7dr/gknaN0LL6iysrLfbbjt7Wp++lF11myQScS7lyw+ZoKK5y+UHS3IYPUAgCNx451q/cuLiu/cLhlXCkVUMGOOSs++THZhYcZed+3atWptbdUvf/UrPfzQQym3k+trF+T8qgnZ6soJhUIq3P8L0dnZOeDrVjc//YRan18luY5iE6YqXF6hZEOdOjdvUufWGslHXVIAkHdsW3LdAz6f69X+6vNqf329ik//gkrPvjgjL3vGGWfo+eefH3Q7rusql9fL9X0YcF1Xx8+bp2XLlumu73639/5nnnlGX7z0Uq1cuVJfvOQSSVJjY6POXrpUW7Zs0T/fdZeqqqr69RrNTz+h1jV/UOGseaq44ApFh4/q3Rbfu0v1T/5K7ZveGNT/AwCQusIZxx3287l1ze8lGZWefUm/2hrIcSVdcj3GIadjBowxg55FYNu2/uGWW3TfffepoaFBkrRhwwZdsWKF7rjjjj5v2NChQ/XqK6+o5t139ev/+i/t2bPnqO277e1qfX6VCmfN04hrburziyZJ0eGjNOKam1Q4a153MgUAZI9t9+vzufX5VXLb2/vZZP+PK+mSjuPhYOTF0evLX/6yKisrde+99+qTTz7RxZdcouVf+Yq+feONh3z8iBEjNLu6WuvWrz9q281PPyq5jiouuOKwVy+0bFsVFyznVAEAZJvr9vPz2VHz04/1u9mBHlf8Luc9A+kQDod180036d6f/UwXX3KJjj/uON199919HrNnzx41NzdLkpqbm7Vu/XpNnTLlqG131mxQbMLUgxLnZ0WHj1ZswtHbAwCkT2zClH5/PnfWvN3vdvtzXEk3egbS4Mtf/rLa29tljNGDDz6oUKjvUIwdO3ZoyZIlOmnBAp31+c/rur/7O82ePfuo7ZpkQuHyin7VEB5aKXnsGtUAkLcsS+Hy/s0ICw+tlEkObOba0Y4rknT+BRfoihUr9NRTT2nS5Ml6/fXXB/QaXpHTAYRWGg+c377pJklSXW3tId+wefPm6ZVXXhlwu1Y4omRDfb8em2ysk3yykiIA+J4xSjbU9euhycY6WeHIgJo/2nFFkn735JMDavNI0nlMHKi86Bm44447tGrVKq1ds0ZJx9EDDzyQtrYLZsxR17b3Fd+764iPi+/dqa5tm9P2ugCAo+vatrnfn88FM+b2u91MHle8KKdhIB0paOXKlfrxT36ixx57THPmzNG3rr9eP/rXf03bQkalZ18m2SHVP/mrw64QZVxX9U8+xGwCADiKtvJR+mTWYu09dp4aR05Se+kwOeFo6g3adj8/n0MqPfvSfjWZ6ePK4eSyZyDnKxDGu7pSnl/51FNP6UuXX64HH3hAF110kaTuwYHTpk/Xv3zve7ryyivTUmPz079R65rfH2adgZ2qf/Ih1hkAgH5IRIfotUtvkwn17bIPxTsV7WhStL1Z0fYmRTuaFdt/O9bWqIKWOoUTnYdt92ifz8WLl/Vr4aFsHVc+y7ZtRWOxjLTdHzkPA8lEIqXLF7/xxhs6e+lS3X777frW9df32XbnnXfq0cce01tvvnnY8zwD1XcFwikKD61UsrGu+9SAbcsuKVOkvH+LGAFAkG2cdqb2Dp884OdFEh0a0tGsIZ3NvX8XdjQrvHubQnW7ZLnuIT6fQ/1egTDbx5UDhcNhhSMDG9OQTjkPA67jKB6P57KEfuu+NsFj3dNTnKSscESF1fNUft7lChcW57o8APCFXc2denZLbVrbDNvSkI5mxWp3aNzG1SrqaFbBjLkqPfvSjF6bIF2i0ajsDISM/sp5GDDGqKvz8F0/XhSJRjOSDAEgCIwx+s2m3WqLp//qgnNHlWrW8CIlfPIls0esoCDYswksy1I4h1dqSoXNQEEASJlrpMrCQQwaPIx5o8s0Z1Sp7z6jw+FwToOA5IELFUlSKBxOadxALnjhTQMAP2ruTOj92jZtrWtX3Env8u3zx5Rp5ogSSZ9+yfTLcSWXly7ukfsK1P3GhUIhOU76u4zSzQtvGgD4heMafdLUofdr27S7pSsjr3HiMUM1fXjfcVt++ZIZCoU88QXTM0e2cCTi+TAQjkQ88aYBgNe1dCW1pbZNW+ra1JnM3EXc5o8pOygISPt7ByIRJTO8NsBg5XIGwYE8EwYsy1IkGvXsoA/Lthk0CABH4BqjHU2der+2TTubMz8w/LhRpb2nBg6lp8f5cAsS5VokGvXMF0zPhAFp/xtn2ykvQpRJUXoFAOCQEo6rzbVtqtnbqvZEdnp4q0eUaPao0iM+xrIsRSMRdXVl5vTEYNge+4LpqTAgdSeleFdXTi/l+FmRSOSw18oGgKDqTDr6695WvbevVXEne5/Z04cV67jRRw4CPSzbViQSyfhSwgPR0xPuJTlfZ+BQjOt6JsmFIxHfTX0EgExq7UqqZm+rNte2ycnyIWRKVZEWjB064J7aZDLpmfEDsVjMc18wPRkGJMl1XSXi8Zz2EBAEAOBTjR0JbdrTog/r25WLT+ZjKwp16vjylE/Z5joQ9PQIeHEdBM+GAal7lapcnTJglUEA6LavtUsb97Tok6bcrRY7bugQnXZshexBjt1yHCcnA9Uty1I0FvPs2DNPhwGpOxA4yWTW5ov2nF/yYnIDgGwxxmhnc6c27mnR3tbczvIaU1qgMyZWKmSn50Dquq4SiUTWZhmEw2GFPL5gnefDQI9svHnhSMQzC0AAQC64xuijhg5t2tOiho7Md6kPLYho6rAiDYmEtPaDuoO2jyyJ6XOTqtIWBHoYY+Q4TkZPG/jpy6VvTojbtq1oNCrXdZVMJtMaCvyQ2gAgk5Ku0da6Nr27p0WtGbiA0IFCljS+vFBTq4pUVdQ9194Yo+JoqM9rVwyJpLVH4EA9SxaHQqG09z5btq1wOCzbtn1zXPFNz8Bnua4rx3HkpPgG2ratkM/eLABIt3jS1V9rW/Xe3taMrhQoSaWxsKYOK9LEiiLFwgd/W96wq1lv72qWJBVHQzpn2nANiWRn7JYxpvu4kkymvNZNaH+48ENPwGf5NgwcqOdNNK7b/fdn/0uWJduyug/8+w/+BAAAQdaecPTe3ha9v69NCTdzhwHbksYOHaKpVUUaUXzkAXRt8aSe2LhbsbCtc6YOV2lBbjqvjTHdf/YfU1xjpM8cV6wDjin58KXSN6cJjqTnQkdi9D8AHFFHwtGGXc3aUtemDGYAxcK2pg8r1pSqon5/uy+KhjV+6BDNGlGSsyAg6dMvjLatoBxV8qJnAABwZEnXqGZvizbtbsloT0BRNKRZI0o0qbJI4RTO9bvGDHr6IAaOMAAAecwYo20NHXpzZ5PaMjgwcGhBRNUjSzS+fAgHcx/Ki9MEAICD7W3t0uufNKquPXPT54YXRTVrZInGlBb4/rx5kBEGACDPtHQl9caOJm1v7MjYa4wpLVD1yBINL45l7DWQPYQBAMgTXUlXG3c36719rRkZHGhJmlBRqFkjSlQ+JJL+F0DOEAYAwOdcY/T+vja9vatZcSf9awWELGlyVZFmDi9RcYzDRj7iXQUAnzLG6JOmTr2xo0nNXem/fks0ZGnasGJNG1actcV/kBuEAQDwofr2uF7/pEl7WrvS3vaQiK2Zw0s0papIkZD/VtPDwBEGAMBH2uOO3trZpK317WlvuzQW1swRJZpYUZiR6wHAuwgDAOADCcfVu3tatGlvq5w0jw4sioZ0/Ogy1ggIMMIAAHiYa4w+qG/XWzub1JFI7+DAiG1p9shSTR9eTE9AwBEGAMCjmjoTevGjBtW2xdPariVpSlWR5o4qVQEDAyHCAAB4jmuM3t3Tord3Nad9vYAxpQWaP6ZMZawTgAMQBgDAQxo7Enrxo/q0LyE8tCCi+ceUaXRpQVrbRX4gDACAB7jGaNOeFm1Ic2/AkLCtuaPLNKmykMGBOCzCAADkWEN7XC9+1KD6jvT1BoQsSzNHFGvWiBLWCsBREQYAIEdcY7Rxd4ve2Z3e3oCJFYU6bnSpiqJ8xKN/+E0BgByo398b0JDG3oDhxVGdMGaoKouiaWsTwUAYAIAsclyjjbub9c7uFqWrM6AkFta8MWUaW1Ygi3EBSAFhAACypK49rpfS2BsQDVmaM6pUU6tYNAiDQxgAgAxzXKMNu5u1KY29AdOGFWnuqDLFwgwOxOARBgAgg2rb4nrxo3o1dabnEsMlsZBOGVehESWxtLQHSIQBAEhJMpmU67qKRg89WM9xjd7e1ax396SvN2DG8GIdN7pUYZveAKQXv1EA0E+//e1vtXz5ck2bNk2RSESxWEynn3664vG+1w7Y19al39fs0aY0BYHSWFhLpw7TCccMJQggIyxjTJpXvgaA/FNTU6OZM2fqpJNO0oIFCzR79my1t7frxhtv1HPPPaczzzxTxhht3NOit3c2pyUEWOruDZg7ukxhBggigzhNAAD98OGHH0qSHn30UTU2NsqyLJWXl+vGG29UV1eX4klX6z+q1ydNnWl5vdJYWAsnlGtYEWMDkHn0NwFAP0yYMEGSNGXKFM2dO1e//OUve7e1dCb0h/f2pCUIWJJmjSjRshkjCALIGsIAAPTDzJkz9fTTT+tHP/qRRo4c2Wfb6580qjXuDPo1ygrCOmfacM0bU8a6AcgqwgAA9NOSJUt0/fXXq7i4uM/97iDbtSRVjyzRedNHqIqlhJEDjBkAgAFK57DroQURLZxQrspCQgByhzAAAAOwo6lDbfHBLyBkSZo9skTVI0s5JYCcIwwAQD+4xuidXc167cO9ct2+Jwa6OjvlOo7sUKhfbZUPiWjh+HJV0BsAj2DMAAAcRWfS0Z+31upbf/c1Xbtoqvbu2N5n+///D1/XTeefoh0fbj5iO7YlzR1VqnOnDycIwFPoGQCAI6hti+v5D+u0Z1+t1v3+Md16662aN2+eZsyYocrKSj3yyCOSpBtuuEGvPvN7Xfz1bx+ynaJoSGdMrGRsADyJMAAAh2CM0ebaNr32SaNcI4VCYdmhkLZv365YLKZNmzb1eXwikVA4cugD/ZjSAp06oYIrDMKzWI4YAD4j6bp6ZXujPqhv73P/71b+VGt/85Ac5+A1BUZNmKRrb79b5cP6rkEwd1SpZo8skWUxSBDeRRgAgAM0dyb1/Id1auhIDKqdWMjWomMrNLq0IE2VAZnDaQIA2O/jxg6t31avhDu470iVhRGdMbFSRVE+YuEP/KYCCDxjjN7a1ayNu1sG3dbUqiKdcMxQ1g6ArxAGAASaa4xe+qjhoPEBAxWyLJ08bqgmVhalqTIgewgDAAIr4bh6/sM67WzuGlQ7JbGwzji2QuVMG4RPEQYABFJHonshobr2wQ0UHFtWoIUTKhQNMW0Q/kUYABA4LV1Jrd5Sq5au1K8xYEk6fnSZZo4oZtogfI8wACBQ6tvjWr2lVp3J1C88XBC2ddqxFRpZwrRB5AfCAIDA2NXcqbUf1A1q6uCwoqhOP7ZShdH+XZQI8APCAIBA+LC+XS9+VK/BLCEwfXix5o8pk81pAeQZwgCAvFezt0Wvf9KU8vPDtqVTxpVrQkVhGqsCvIMwACBvGWP05s4mbdrTmnIbQyK2zppUxbRB5DXCAIC8lI7FhEpjYZ01uUrFMT4qkd/4DQeQd9KxmFBVUVSfm1SpWJiBgsh/hAEAeaUz4ei5QS4mNKa0QKcdW6EICwkhIAgDAPJGOhYTmlRZqJPHlTNjAIFCGACQF+rb43puS606BrGYUPXIEh03qpQVBRE4hAEAvrerpVNrtw5uMaETjxmq6cOL01gV4B+EAQC+9nFjh57/sC7lxYRsS1o0oULjy1lDAMFFGADgWzuaOgcVBCK2pcWTKrnGAAKPMADAl3a3dGrtB7UpB4EhEVufm1SlChYTAggDAPxnb2uX/ry1Tk6KQYDFhIC+2BMA+Erd/lkDyRS7BCoLI/rc5CoVsJgQ0IswAMA3GjoSenZzbcqzBkaXFuh0FhMCDkIYAOALTZ0JPbt5n+JOausITKoo1MnjWUwIOBTCAADPa+lK6tnNtepMcUGh6hElOm40iwkBh0MYAOBpbfGkntm8T+0JJ6Xnzx9TppkjStJcFZBfOHEGwLM6Eo6e3VyrtnhqQeD40aUEAaAfCAMAPKkr6ejZzfvUnOJFh2aPLFH1yNI0VwXkJ8IAAM+JJ109u7lWjZ2pBYGZw4s1dxRBAOgvwgAAT0k4rp7bWqv6jkRKz59aVaR5Y8oYLAgMAGEAgGckXaM/b63TvrZ4Ss+fVFGok8YOJQgAA0QYAOAJjmu09oNa7WntSun548uH6OTx5QQBIAWEAQA55xqjFz6s087m1ILAMWUFWjShggWFgBQRBgDklGuM1m+r18dNnSk9f1RJTKcfW0kQAAaBMAAgZ4wxemV7g7Y1dKT0/BHFUS2eVKmQTRAABoMwACBnNu5u0Za69pSeW1UU1ZmTqhS2+RgDBou9CEBObG/o0Fu7mlN6bsWQiM6aVMXVB4E0YU8CkHX17XGt+6g+peeWFYR11pQqRcN8fAHpwt4EIKs6Eo7+vLVOjmsG/NySWFifnzJMBeFQBioDgoswACBrHNdozda6lK5AWBQNacmUKhVGCAJAuhEGAGSFMUYvfdSg2vaBry44JGJryZRhKopy1XUgEwgDALJi454Wfdgw8JkDsXB3ECiJEQSATCEMAMi47Y0demvnwGcOhGxLZ02uUllBJANVAehBGACQUfXtca3fltrMgUXjK1RZGE1zRQA+izAAIGN6Zg4kU5g5cNyoUo0rH5KBqgB8FmEAQEY4rtGaD1KbOTChfIiqR5ZkoCoAh0IYAJB2xhi9tL1BtW0DnzlQWRjRKeMruBQxkEWEAQBpt2lPiz6sH/jMgcJISIsnVSnMhYeArCIMAEirjxs79GYqMwcsS4snVbKoEJADhAEAaVPfHte6FGcOnDqhnJkDQI4QBgCkRUfC0ZoPUps5MHdUqcaXF2agKgD9QRgAMGiOa7T2gzq1xVObOTCbmQNAThEGAAyKMUYvb2/QPmYOAL5FGAAwKDV7W/VBqjMHJjJzAPACwgCAlO1t7dIbO5oG/LzemQNRZg4AXkAYAJCSrqSjFz6s18CHCzJzAPAawgCAATPGaP22hpSWGp7DzAHAcwgDAAbs3b2t2tHcOeDnjR86RHOYOQB4DmEAwIDsa+3SmymME6gsjGjhhHJmDgAeRBgA0G9dSVfPpzBOYEjE3j9zgI8cwIvYMwH0izFGL35UP+BxArYlnTGRmQOAlxEGAPRLzd5WfdI08HEC88aUaVhRLAMVAUiXcK4LAJB9xknKqd8j09UhOcnuO0NhWbEhClWMkBXq+9Gwry219QTGlhVo+rDidJQMIIMIA0AAuG3Nin+wUU7tDiX3fCy3cZ9kDnPm37JkDx2m8IixClWNUXRitXY0mQGPEyiKhlhqGPAJy5jDfSIA8DPjukp+vFmdm15Wcvt7kpFkWZJx+9eAZXcHBkuKjJ+uvZNO1WsdQxR3jv6RYVvS0qnDVVXEwkKAHxAGgDxjjFH8/TfU8fJTMu3N+w/q/QwAh7O/jY6yUXp75vmqdY98kD/hmDLNGM56AoBfEAaAPOI016t9zRNK7tiSsddwZemDGWfrr6VTDnnq4JiyAi2eWMnpAcBHCANAHjDGqOudF9Xx8p8k1x18T8DRWLbqS0frrenL1H7A0KOiaEjnTR+hWJiJSoCfEAYAnzOuo/a1v1H8vdez/trxUEyb5l6iHZFKWZKWThvGNELAhwgDgI8Zx1Hb079SYtu7uatB0s7pn1d45kmaPqI0Z3UASB1hAPAp47pqe/ZhJbZulFK6kHB6RSbNUdHnvyyLJYcB32GvBXyq8801Smx9R14IApKU2LpBnW+uzXUZAFJAGAB8KLlvhzpfeybXZRyk87Wnldy3I9dlABggwgDgMyaZUNszD0vy4tQ9S23PPCyTTOS6EAADQBgAfKbjtWflNtVlfvpgKowrt6lOHa8/m+tKAAwAYQDwEberQ13vrJdXxgkcmlHXhvUyXQO/wiGA3CAMAD4Sf+/1T68y6GVOUl05WPcAQGoIA4BPGNdV54b1uS6j3zo3rJNxPXgqA8BBCAOATyQ/fl+mtTFrr9ceT2jObT/T/3z8uZSeb1oblfz4/TRXBSATCAOATyQ+2SJlcUGfH/3pJc2fMCr1BmxbiQxeMAlA+hAGAJ9I7tnefRGiLNi6t17v76nTklmTUm/EdZXcvT19RQHIGMIA4APGdeTU7hxUG65rtOCO+3T7b9b0uX/1ux9oxA0/1H+/8V7vff/riT/rf114xqBeT5Kc2p2MGwB8gDAA+IDbsG/Qswhs29KNS0/W/S+8qcb27ml/Gz/Zq6t//lvddsHpumjedEnSH9/erEnDKzR5RMWg65aTlNuwd/DtAMio8NEfAiDXnOb6tLRz2Ymz9IM/rNf//fPrWrFwji6/91FdftIs3bBkQe9jXv9wp554vUa/feM9tXUllHAclQyJ6dZzT02t9pZ6hSpHpqV+AJnBVQsBH4hv2aC2Zx5KS1sPvPCW/veTazWqrETjKsv0n393sUKHGZj40EvvqGbnPv3vL34u5dcrWrJc0clzUn4+gMzjNAHgB2k8737piTPVEU/KyOi+a84/bBBIF+M6GW0fwOBxmgDwg1AobU39j//qvtphfWvHUYPA8lNmD/r1rHBk0G0AyCx6BgAfsKKxtLRz1++e19Mbt+rpf/gbJV1Xv3xxQ1raPRIrkp7aAWQOYQDwgVDlIBb/2e8/1r+te599TQ9944uqPma4rjvzBP3kmVeUcDLbjc/gQcD7CAOAD9iFJbKGFKf8/Gc2bdWtv35G/+eqZTrx2DGSpGsXz1dLR5d+/cqmdJV5EKuwRHZhScbaB5AehAHAJ8IjxkmWNeDnvbV9t675+W91+0WLdf7x03rvLx0S07WL5+vHT78sJxMLA1lWd80API+phYBPdLzxZ3W++rTkl13WsjTkpKUqmLc415UAOAp6BgCfiE6a7Z8gIEnGKDJp8LMRAGQeYQDwiVBZlcLjpkmWD3Zby1Z43DSFyipzXQmAfvDBpwqAHgWzT5WMDy78Y1wVzElt+WIA2UcYAHwkPHay7NKKlAYSZo1lyS6tUPiYybmuBEA/EQYAH7EsW4WLLvT22AFjVHjahbL8cDoDgCTCAOA7kfHTFJ15kjd7ByxL0ZkLFBk37eiPBeAZhAHAhwoXnie7uMxbgcCyZBcPVeHC83JdCYABIgwAPmRFYir6/FckeSsMFH3+K7Ii0VxXAmCACAOAT4VHjlfR0hXe6B2wLBUtXaHwSFYcBPyIMAD4WPTYmSo6+wrJtnMTCixLsm0VLV2h6ISZ2X99AGnBcsRAHkjs/ECtf3xQSiaytw6BZUvhiIrPvVKR0ROz85oAMoIwAOQJt61Z7c//RoltNVl5vciEmSo8/SLZRaVZeT0AmUMYAPKIMUaJDzaq/fnfyHR1pH89AsuSFRuiwtMv7r5WAoC8QBgA8pDb2a6ujS+pa9PLMu0t3V36qZ4+2P9cq7BEsVknKzb7FNmxwvQWDCCnCANAHjOuq+TH76tz08tKfvRXSft39yOFgz7bLIXHT1PBrJMVHjtVls2YYyAfEQaAgDCJuJy6XUrW7pSzb4ecfTvkdrVLyWT3A8Jh2bFChYaNUWjYGIWrRitUOYp1A4AAIAwAABBw9PkBABBwhAEAAAKOMAAAQMARBgAACDjCAAAAAUcYAAAg4AgDAAAEHGEAAICAIwwAABBwhAEAAAKOMAAAQMARBgAACDjCAAAAAUcYAAAg4AgDAAAEXDjXBRwoevw1ssNRWXZIlh1SKPLpbcu2P90WCskOR2X3bgsdtM2yQ7JtS5ZtKRSyZX3mtm1bskNW72OOuM2yFArbCtmWQral6P7b4d5/hz7dFvr0ceEDHhs61G3Lkm1ZCllSJGT33g6HbIUsdf/bthSxrUPc7t4ese3e2yHLkmVJtiVZlva3L1mSQrYlW+r+v9jqvW1bUsg68HZ3G5YxknFluUmpz223+497+G2WcSXH+fS2m5RcR8Z1pWRcxnEk1+2+L5mQcZ3u24mE1HO757E9j0vEP32O68hNJGUcV8Z15caTcp3u5xjHlZtIynU+vW3233YSSZkDHufEkwfcdmRcI9cx+/+9//mu6d7mGBnHyHVcOQl3f5tGTsLZ/5xPn+caI8cYxV0jx+gztz/77+7brrpvO0b7t316+/+YbTndL9OF/Zv9m/3bu/s3PQMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAASdyVOdnZ3mO9/5juns7Mx1KQfxcm3GUN9geLm2fOLln7OXazOG+gbDy7UNlmWMMbkOJJnQ3NyssrIyNTU1qbS0NNfl9OHl2iTqGwwv15ZPvPxz9nJtEvUNhpdrGyxOEwAAEHCEAQAAAo4wAABAwOVtGIjFYvrOd76jWCyW61IO4uXaJOobDC/Xlk+8/HP2cm0S9Q2Gl2sbrLwdQAgAAPonb3sGAABA/xAGAAAIOMIAAAABl3dh4JZbbtFpp52mK664QvF4vM+2jo4OLVu2TGeccYaWLFmi+vp6T9XX43vf+55OOOGEnNeUTCZ11VVX6bTTTtPf//3fZ62e/tbXI9s/rwMdrjYv/K7lI/bv9NXE/n10Qdq/8yoMvPnmm9q9e7deeOEFzZw5U4899lif7X/6059UXV2ttWvX6ktf+pL+8z//01P1SVJLS4s2btzoiZp+97vf6ZhjjtELL7yg9vZ2vfjii1mrqz/1Sdn/efW3tlz/ruUj9u/01sT+nXptuf5dy4S8CgMvvfSSzj77bEnSOeecc9Av95QpU9Te3i5Jamxs1LBhwzxVnyT9+Mc/1vXXX++JmvpTby7rk7L/8zrQkWrL9e9aPmL/Tm9N7N9HFrT9O5zrAtKpsbFRo0ePliSVlZUd1HUzadIkbdy4UdXV1bIsS6+88oqn6mtqatI777yj2267zRM1NTY29q6/fah6c11fLn5e/a0t179r+Yj9O701sX+nXluuf9cywZc9A7t379aiRYsO+mOMUXNzs6TuN7KioqLP8x588EEtXrxYGzdu1B133KE777zTU/Xdc889+ta3vpWRmg6nvLz8sDUdaZsX6svFz+tAR6otW79r+Yj9O33Yv1MXtP3bl2Fg5MiRWrdu3UF/zj33XD399NOSpKeeekqnnnrqQc/teUOHDh2qxsZGT9W3ZcsW3XXXXTrnnHO0efNm/cu//EtG6jvQySeffNiajrQtW45UQy5+Xv2tTcrO71o+Yv9OH/bvzNQm5eH+nburJ2fGzTffbBYtWmSWL19uurq6jDHGfP3rXzfGGNPU1GTOPfdcc8YZZ5hTTz3V/PWvf/VUfQeaP39+zmrqqSeRSJivfvWrZtGiReaGG27IWj39re9A2fx5HehwtXnhdy0fsX8Pvib27/4L0v7NcsQAAAScL08TAACA9CEMAAAQcIQBAAACjjAAAEDAEQYC4IEHHtDQoUPT0ta2bdtkWZbC4bB27NjRZ9uuXbsUDodlWZa2bdvWZ9vjjz+uxYsXq6ysTMXFxZozZ47uvPPO3oU80lkjEDRXXXWVLMvSddddd9C2b37zm7IsS1dddVXvfbt379YNN9ygiRMnKhaLaezYsTr//PO1evXq3sdMmDBB99xzTxaqhxcQBpCS0aNH6z/+4z/63Pfggw9qzJgxBz32n/7pn3T55ZfrxBNP1J/+9Cdt3LhRd999t95+++28WNMb8IKxY8fqkUceUUdHR+99nZ2devjhhzVu3Lje+7Zt26b58+frueee0w9+8AO98847WrVqlc4888ycLf2L3CMM+MCqVau0aNEiDR06VJWVlVq2bJm2bt0qSVqzZo0sy+qz6MVbb73V++18zZo1uvrqq9XU1CTLsmRZlm6//XZJUkNDg7761a+qvLxchYWF+sIXvqDNmzf3q6Yrr7xSK1eu7HPfAw88oCuvvLLPfa+++qr++Z//WXfffbd++MMfauHChZowYYKWLFmixx9//KDHA0jNvHnzNG7cOD3xxBO99z3xxBMaO3asjj/++N77enoKXn31VV166aWaOnWqZs2apZtuukkvv/xyLkqHBxAGfKCtrU033XSTXnvtNa1evVq2beviiy+W67pHfe7ChQt1zz33qLS0VLt27dKuXbt0yy23SOruWnz99df15JNP6qWXXpIxRueee64SicRR273gggvU0NCgdevWSZLWrVun+vp6nX/++X0e96tf/UrFxcX65je/ech2ODUApM/VV1/dJ6Tff//9uuaaa3r/XV9fr1WrVun6669XUVHRQc9nfwyuvLpQUb764he/2Offv/jFLzR8+HC9++67R31uNBpVWVmZLMvSyJEje+/fvHmznnzySa1fv14LFy6U1H3gHjt2rP77v/9bl1122RHbjUQiWrFihe6//34tWrRI999/v1asWKFIJNLncZs3b9bEiRMPuh9A+v3N3/yN/vEf/7F3bM/69ev1yCOPaM2aNZK6l/g1xmj69Om5LRSeQ8+AD2zdulXLly/XxIkTVVpaqmOPPVaStH379pTbrKmpUTgc1oIFC3rvq6ys1LRp01RTUyNJ+sIXvqDi4mIVFxdr1qxZB7Xxt3/7t3r00Ue1e/duPfroo32+gfQwxsiyrJTrBNB/VVVVOu+88/Tggw9q5cqVOu+881RVVdW7vWfBWfZJfBY9Az5w/vnna+zYsbrvvvs0evRoua6r6upqxeNxFRcXS/p0J5fUr27+w61CfeDB++c//3nvYKRDfbOvrq7W9OnT9ZWvfEUzZsxQdXW13nrrrT6PmTp1qtatW6dEIkHvAJAF11xzTe/V/n7605/22TZlyhRZlqWamhpddNFFOagOXkXPgMfV1dWppqZGt912m8466yzNmDFDDQ0NvduHDRsmqXtaX4/PHpCj0agcx+lz38yZM5VMJvtch7uurk7vv/++ZsyYIUkaM2aMJk+erMmTJ2v8+PGHrO+aa67RmjVrDtkrIEnLly9Xa2ur7r333kNuz4urfQEecs455ygejysej2vp0qV9tlVUVGjp0qX66U9/qra2toOey/4YXIQBjysvL1dlZaX+/d//XVu2bNFzzz2nm266qXf75MmTNXbsWN1+++16//339Yc//EF33313nzYmTJig1tZWrV69WrW1tWpvb9eUKVN04YUX6tprr9W6dev09ttva8WKFRozZowuvPDCftd37bXXat++ffra1752yO0LFizQrbfeqptvvlm33nqrXnrpJX300UdavXq1LrvsMj344IOp/WAAHFIoFFJNTY1qamoUCoUO2n7vvffKcRyddNJJevzxx7V582bV1NToJz/5iU455ZQcVAwvIAx4nG3beuSRR/SXv/xF1dXV+va3v60f/vCHvdsjkYgefvhhvffee5o7d66+//3v67vf/W6fNhYuXKjrrrtOl19+uYYNG6Yf/OAHkqSVK1dq/vz5WrZsmU455RQZY/THP/5xQN354XBYVVVVCocPf8bp+9//vh566CG98sorWrp0ae80pjlz5jC1EMiA0tJSlZaWHnLbscceqzfeeENnnnmmbr75ZlVXV2vJkiVavXq1fvazn2W5UngFlzAGACDg6BkAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgID7f9cBN5rF5EHjAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tau_max = 2\n",
    "pc_alpha = 0.01\n",
    "\n",
    "pcmci = PCMCI(dataframe=dataframe,\n",
    "            cond_ind_test=ParCorr(),\n",
    "            verbosity=0,\n",
    "            )\n",
    "results_pcmciplus = pcmci.run_pcmciplus(tau_max=tau_max, pc_alpha=pc_alpha)\n",
    "tp.plot_graph(\n",
    "    graph = results_pcmciplus['graph'],\n",
    "    val_matrix= results_pcmciplus['val_matrix'],\n",
    "    var_names=dataframe.var_names,\n",
    "    ); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here PCMCI+ misses the link $X^2\\to X^3$ and fails to orient the link $X^3\\to X^1$. \n",
    "\n",
    "### Bagged-PCMCI+\n",
    "\n",
    "From our numerical experiments, we found that Bagged-PCMCI+ improves upon precision of adjacencies and recall of orientations. However, one cannot directly compare results from Bagged-PCMCI+ versus PCMCI+ for the same ``pc_alpha`` (see precision-recall curves in the paper). Here we choose a higher ``pc_alpha`` for Bagged-PCMCI+. Further below we illustrate how to use model-selection to choose ``pc_alpha``.\n",
    "\n",
    "Another parameter of Bagged-PCMCI+ is the block-length of the bootstrap (by default 1). It can optionally be used to better deal with autocorrelation, but its effect was not yet evaluated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAHBCAYAAAD0E7h1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9jElEQVR4nO3deXwc9Z3n/3dVdbduybIkH/JtywZscxsMtgmQACYEcpAEEsKEY0IWAnlkgGwy2WEXyOTOMstkJ2R+k4RjJgHPcGR+yRAgxMSAOcxhsDG28X1btnVZt7q76rt/tCUsbOHuVre6qvv1fDz8sNxdVf2xWqV697e+h2WMMQIAAAXLznUBAAAgtwgDAAAUOMIAAAAFjjAAAECBIwwAAFDgCAMAABQ4wgAAAAWOMAAAQIEjDAAAUOAIAwAAFDjCAAAABY4wAABAgSMMAABQ4AgDAAAUOMIAAAAFjjAAAECBIwwAAFDgCAMAABQ4wgAAAAWOMAAAQIEjDAAAUOAIAwAAFDjCAAAABY4wAABAgQvlugAAI8sYk/jjeTKJBxJPWJYsSZZty7IsWZaVwyoBjCTCAJDnjDHyPE+e68rzPJn+i/8xWJYl27ZlO47sQwEBQH6yTLK/GQAEijFGbjwu13WTDgBDsSxLjuPICYUIBUAeIgwAecYYo3g8Ljcez8rxnVBIIUIBkFcIA0AecV1XsWh0RF4rHInIcZwReS0A2UUYAPKAMUaxWEye647o69qOo3A4TCsBEHCEASDgjDGK9vUNu19AuizLUqSoiEAABBjzDAABlusg4JcaAAwPYQAIKD9dhP1UC4DUEQaAgIrFYr66+Pb3WwAQPIQBIIBc1x3xzoLJ8FxXrg/rAvDhCANAwBhjRmz4YDpi0aivWiwAHBthAAgYPweBfkGoEcD7CANAgPSvM+B3qayBACD3CANAgMSzNMVwNgSpVqDQEQaAgOhfeCgo3Hic1gEgIAgDQECM1OiBnbt26aLFi3XqaafpjDPP1ONPPJH2sfw44gHAkZiOGAiIWDQ6IsP29u7dq/379+vkk0/W/v37dfaCBVq9apXKyspSPpbjOApHIlmoEkAmhXJdAIDkjFTHwfHjx2v8+PGSpDFjxmj06NFqaW1NKwwEobMjAG4TAIFgjBn2/XfP83TyKafo7+64Y9Djzz77rCqrqo56O+DNN9+U53maNHFiWq+ZiboBZB9hAAiATFxQbdvWf//mN/XLX/5Sra2tkqTVq1frS1dfrbvvvlufvfzyQds3NzfrKzfcoH/6p38a1usSBgD/o88AEACu62ZkIp94PK4TTzpJV3/pS7rmmmt07nnn6bJLL9W99947aLu+vj594tJLdf111+mqq64a1muGIxE5jjOsYwDILsIAEABuPJ6xRYB+9atf6c677lJ9fb2mTJ6sf//3fx90sTbG6Jprr9WsmTN1xwduKaQjHA7LCdE9CfAzwgAQAJkMA52dnZo0ebJmzJih55ctO6Jj4Esvv6wLL7xQJ86dO/DYr3/9a8097N+pIAwA/scZCgSBZWXsULfedpskqbmp6ajN9wsXLFB3V1fGXi+TtQPIDjoQAgXk7rvv1tNPP63nly1T3HX14IMP5rokAD5AGAACwLaHf6o+8MAD+sef/UyPPfaYTjrpJN1y8836h//zfzJ2+2EomagdQHZxlgIBYFnWsJrbn3nmGf3Nrbfq/l//WvPPPFOSdNNNN6mjo0MPP/xwpso8kmUlagfga4QBICDS/YS9cuVKfenqq/X9739fn/70pwcer6ys1E033qj/fc89WZvmmFYBIBgYTQAERDweVzzLTfqZFgqHFWIkAeB7xHYgIII4cU8QawYKEWEACAjLsgLV7G7bNv0FgIAIzm8WAAqFw7kuIWlBqhUodIQBIECC8mk7aK0YQKHjbAUCJhyJ5LqEYwpCjQDeRxgAAsa2bV/30A+FQrQKAAHDGQsEkBMK+fJ2gWVZLEoEBBBhAAggy7J82RQfjkR8GVIAfDjCABBQtm0rUlSU6zIGRIqKuD0ABBRnLhBgtm0rbOJSZ0tO6yAIAMHGzT0gwIznymx5U05Xm8yYGfLGTpeskbso99+uIAgAwcbaBECAudtWy9u35f0HykbJnXSSFCnJ+ms7oZBCPu3ICCA1hAEgoLyWPXI3vnbkE05I3oQ5MqPGZeV1aQ0A8g9hAAgg09ul+Jq/SG586I1qJsmrP14mQ7cNrEPzGwRlFkQAyaPPABAwxnPlbnr9w4OAJDXvVKhitKy6KXJdV67rSqlmf8uS4zhyHIeWACCPEQaAgPF2vivT1XbM7axRY2WPmTqwTkA4HJYxRsbz5B362xij/nhgKXELwLJt2Yf+pgUAKAyEASBAvJY98hq3HHvDSImc6acdcTG3LEuW4zCmGMAg/E4AAsL0dsndsjKJLS05DfNkhf0zIREAfyMMAAFgPC+5fgKS7EmzZVfUjEBVAPIFYQAIAG/7O8n3ExjfkP2CAOQVwgDgc27jFnn7tx57wyH6CQDAsRAGAB/zDu6Xt/2dJLaknwCA9BEGAJ8yPR2HZhg89twA9BMAMByEAcCHTDyq+HuvJtVhkH4CAIaLMAD4jPG8RItAX9exNy4qpZ8AgGEjDAA+YoyRt321THvTsTe2QwrNOot+AgCGjTAA+Ii3b4u8/duS2tZpmCertDK7BQEoCIQBwCe8tn1JjhyQ7MlzZVdnZ4liAIWHMAD4gOluT8wwmASrborscTOyXBGAQkIYAHLMxPoU35DkyIGKGjlTT6bDIICMIgwAOfT+yIHuY29cVCpn5pmybE5bAJnFbxUgR4wxcre9LdPRfOyNHUYOAMgewgCQI17jZpkDO5La1mk4g5EDALKGMADkgNfaKG/HmqS2taecKHvU2CxXBKCQEQaAEZYYOfBGUtvaY6bKHjs9yxUBKHSEAWAEDYwc8JIYOVBZK3vKSYwcAJB1hAFghBjPlbthRZIjB8oYOQBgxPCbBhgBxhi5W9+W6Ww59sZOSKHjzpIVimS/MAAQYQDIOmOMvJ1rZZp2JrG1lWgRKKnIel0A0I8wAGSZt/s9eXs3JrWtPfVE2VVjslwRAAxGGACyyN2zQd7u9Ulta4+dJoeRAwBygDAAZInbuFnezrVJbWtV1smecmKWKwKAoyMMAFng7tua9HLEKi6XM/MMWRanI4Dc4LcPkGHege3ytq1KbuNQJLHmACMHAOQQYQDIIK9pp9wtbyW3sRNW6PiFskrKs1sUABwDYQDIEK9lt9zNK5Pb2AnJOX6BrLKq7BYFAEkgDAAZ4LU2HlpvwBx7Y9uRc9wC2eXVWa8LAJJBGACGyWvbL3fja5JJNgicLbtidPYLA4AkEQaAYfDaD8jduEIy3rE3tmw5s+bLrqzNfmEAkALCAJAmr6NZ7nuvSp577I2txDTDzC4IwI8IA0AavM5Wue+9klwQkCWn4QzZ1eOyXhcApIMwAKTIdLXJXf+y5MaT2t5pOF326PosVwUA6SMMACkw3e2Kr39ZcmNJbe9MP012zcQsVwUAw0MYAJJkejoVX/+SFI8mtb0z7RTZdZOzXBUADB9hAEiC6e1SfN1yKdaX1Pb2lBNlj5ma3aIAIENCuS4A8DvT05G4NRDrTWp7e9IcOeNmZLkqAMgcwgDwIbyOZrkbXpXiyfURsCceL6d+ZparAoDMIgwAQ/Ba9iSmGE5mQiFJdv0s2fXHZbkqAMg8wgBwFG7jFnnbVye9vT1uhuyJJ8iyrCxWBQDZQRgADmOMkbdzrby9G5Pexx4zTfbkuQQBAIFFGAAOMZ4nd+tbMk07k97Hqpsse+pJBAEAgUYYACSZeEzuxtdk2g8kvY89ZqrsqScTBAAEHmEABc9EexR/7xWpuz3pfeyJJ8iun0UQAJAXCAMoaANzCER7ktvBsg7NLDglu4UBwAgiDKBgDSxBnOQ6A7KdxDLEo8ZmtzAAGGGEARSkVOcQULhIoePOllU2Kqt1AUAuEAZQcFKdQ0DF5YkgUFyWvaIAIIcIAygY6cwhYJWPljNrvqxwURYrA4DcIgygIBjPk7tlpUzzrqT3sarHyZkxT5bDaQIgv/FbDnkvvTkEpjGZEICCQRhAXmMOAQA4NsIA8pbX3pQYMRDrTW4Hy5Iz7VTZdZOzWxgA+AxhAHnHGCNv7yZ5O9dKMsntZIcOzSEwJqu1AYAfEQaQV0w8KnfzSpm2xuR3Yg4BAAWOMIC84XW2yt30utTXnfxOxeUKHX+2rCLmEABQuAgDCDxjjLz9W+VtX5P8jILqn0PgLFnhSBarAwD/Iwwg0Iwbk7v1bZnm3SntZ1WPl9MwT5btZKkyAAgOwgACa+u61Xrz6d+pr6tD0+rH6ORZ0xQJh4+5nz1uhuzJcxk6CACHEAYQOE8++aS+/c3b9O76DYMenzimRj/6+rW64sJFR9/RCcmZfprs0fUjUCUABIed6wKAVBjP1ReuvEJ14+r1u9/9To2NjWpra9OKFSs0/5zzdPX/vEf/uezVI3csrVJo7vkEAQA4CsIAAsP0dir+7vPq6u7RpZdeqm3btunqq6/WRRddpOXLl+uxxx7T4sWL9f37/2PQfvaYqQrN+QirDgLAECxjTJKzsgC54zXvlrvlLcmL68Qrb9GGHXvk2LYWn32aYvG4/vTqW1qxYoU2b96sq666Sgf+/BtVVVbKmXaK7NpJuS4fAHyNPgPwNeN58naskbdvy8Bjf/zZXXru9VW64MxTNGFMjf740hv606tvKR6Pq7c3MfVwqKwq0RpQWpmjygEgOGgZgG+Zvm65G1+X6WodcpsNO3brvK/+D50y70w988wzOv/889XX3qpX3ljJ0sMAkCR+W8KXvNZGuVvelOKxIbfZsrtRF99yp+rGT9CSJUv0wx/+UC+88IKeffZZggAApIDfmPAV48bl7Vonr3Hzh263bc9+Lb75f6l0VI2WLl2qJUuW6I477tDdd9+tCy64YISqBYD8QBiAb3gHD8jd+tYx1xbYua9Ji2/5XwqVVeq5557Tk08+qa9//etauHChPv7xjysWiymcxORDAIAE+gwg50w8KnfHGpkDO5La/iM3/K32dUb1/PPPa+nSpbr++ut1+I/x3LlztWrVKtk2I2cBIBm0DCCnvJY9cretkmJ9SW3f3tWtV995T0uWLNGkSZN05plnatWqVQPPv/nmm7ruuuu0efNmzZw5M1tlA0BeIQwgJ0y0V+621TKte1LarzgSVkVZqe699149+uijRzzf1NSkUCik6urqTJUKAHmPdlSMqMRyw9sVX7005SAgScXTTtJD//Yb1dXVqaur64g/5eXluv/++1VbW5uF6gEgP9FnACPG9HYllhtuP5D6zqEiOQ3zZFfVZb4wAChw3CZA1hlj5DVulrdrneS5Ke9vjZ4gZ+qJssLFWagOAEAYQFaZ7na5W9760FkEhxQuljPtZNnV4zNfGABgAGEAWWE8V96eDfL2bJDSuBNlj5kqe9IcWSHmCwCAbCMMIOO8jpbE5EE9HanvXFQmZ/qpsivpAAgAI4UwgIxJTCW8Vl7jlmNvfARL9vgG2ROPl2U7Ga8NADA0wgAywmvbn2gNiPakvnNppULTT5NVNirjdQEAjo0wgGExvZ1yd66VaUl9zgBZtuwJx8se3yCLqYMBIGcIA0iLifXJ271e3v5taXUQtCpq5Ew7RVZJReaLAwCkhDCAlBg3Lq9xk7w9myQvnvoB7JDsyXNkj5kqy7IyXyAAIGWEASTFGE/e/u3ydq9PelGhD7JGjZUz9WRZRaUZrg4AMByEAXwoY4xM6165O9dKvZ3pHSQUkTPlRFk1E2kNAAAfIgxgSF5Hs7wd78p0tqR9DKtmYiIIhIsyWBkAIJMIAziC6elIjBBo3Zv+QSIlcqaeLLt6XOYKAwBkBWEAA0y099AIge2S0lzM0nZkj5she/xMphIGgIAgDEAmHpO3d6O8xs1prSrYz6qbImfi8bIiJRmsDgCQbYSBAmY8T97+rfJ2vyfFo2kfxxo1Ts6k2bJKKzNYHQBgpBAGCpAxRqZlT2KEQF9X2sexyqoTcwawqBAABBphoIAYY2Ta9snbvV6mqy39AxWXJVoCqusZKggAeYAwUACM58o07ZS7d7PUm8aywv1CRbInHi+7bgprCQBAHiEM5DET60v0CWjcKsXTmzVQUmKEwPiGxIJCDiMEACDfEAbykOntlLd3s7ymHcMaHSBZssdMlT3hOFmR4ozVBwDwF8JAnjDGyHS2yNu7Uaa1cdjHs6rrE/0CSsozUB0AwM8IAwFnjCfTslfe3k0yXa3DPp5VUSN70hzZFaMzUB0AIAgIAwFl3Li8A9sTEwX1dQ//gMUVcibPljVqHCMEAKDA5FUYMMbIeJ6MSUylayRZkmRZsg77E2Qm2iOvcYu8/dskNzb8A4aLE7MG1k2WZTFCAAD6GWMG/siY968pUuJ6YtuBv6b0C3QYMMbIjcfleZ48z0tqH8u2ZVuWbMeRHaA30nQflLt3k0zzLsmkuW7A4Uoq5YxvSCwrzDBBAJAxJnE9cV15h32wPBbbtmXbtpxQKDDXlA+yTLL/W5/of7P6Q8BwWJYlx3F8+wYaY2QO7k/0B2g/kJFjWlVjZI9rkFVV58v/MwCMtP4Plq7rJh0AhtIfCoL0YVMKWBhwXVexaPpz6H8Y23EUDoez9uYZY5I+tonH5DXvkrdvq9TTPvwXtyxZNRMTLQGlVcM/HgDkAWOMYrGYPHc4Q7CHFo5E5DhOVo6daYEIA8YYxaLRYbcEJCMbb57XvFsmHpUzdtqQ2xhjZDqa5R3YLtOyZ5jzAxzihBPzBIybzkqCAHCYbH64PJxt2wpHIr5vJfB9nwE3HlcsloGOckmKRaNyM/jmea2Ncje/Iat0lHSUMGCiPfIO7JB3YMewFg0apKhU9rgZiWmDHd+/xQAwYkbyw6UkeZ6nvt5ehcNhOSH//j72bWXGGMXjcbnx+Ii/tud5ivb1KVJUNKxA4B08IHfja4leqF2tMj0dskoqEiMe2hoTrQBt+zJWt1VWnZgyeDQLCAHABxljFO3rG3a/gHTEYjF5xijk0z5qvrxNYIxRPBaTm6X7OKkoKipKq7e919Eid/1Lg5r7rbopskJheQd2Dm+tgA+wqscnQkD5aF/+kAFArhnPU19f5n7vpstxHIWy2D8tXb4MA7FYLCctAkMpKi5O6Y0zXW2Kr3spM/MADMV2ZNdOlj1+hqxipgwGgKEYY9TX25vrMgY4oZDCYX8t+ua72wSu6/oqCEhSNBpVJMk+BKanQ/H1L2cvCISLZI+dLnvMNFnhSHZeAwDyhDFG0RHoKJgKNx5PDEH00UgDX4WB/o4dfmM8T67rKnSMzh+mtyvRIhDPwv+hpELOuBmyaifJsv3zAwQAfua6rswIdRZMRSwalZ1iq3M2+SYM+DUI9IvHYgOzTB2N6e1MBIFYBpuibEfW6PrEqICKGt/80ABAEHiep/gIjkZLVSwa9c2wQ9+Egf7pH/0sFo0edYSB6e1UfO3yjAUBq2yUrLopsmsmygr5674SAASB3z9gShqY+tgPQw5zX4HeH0bod/1TIR9+n8f0dGSmRcAJy66dlGgFKGOWQAAYjlTWFsileDwu23Fy3jrgizAQlDdNSnT86A8DmQgCVmWd7DFTZFWPpy8AAGSI3zqiD+VoHzJzwRdhIChvmqSBFRKtvi7F1y2XYumNW7WqxsiZdrKsorIMVwgAhS2VlWz94PAPmbmS87VrzQi9aR0dHVq4aJHmz5+veWecofvvvz/tY7ldB4cVBCTJxPoIAgAwhLan/0Pd77yW1kiAkZyw7oorr9T4+np98aqr0j6G53k5H/GQ80mHRmrtAdd11dfXp9LSUnV3d+v0efO0/MUXVVNTk/KxLEkRJ9FxUL1dMr2dMof+Vm+XZJJ7U0Mnns8qggBwFE2/+Uf1blyjUN14VZ57qUrmzEt6Nti+3t4Ru/X8/PPPq7OzU7/57W/1yMMPp32cXK9dkPPbBCPVlOM4jkpLSyVJvb29w1q32khSuFh2pESqrBv8nDFStEemt0vq65aJdsv09UiH/a1Dr+s17ZQzmTAAAEOJH9irlsd+qdDz/5VUKDDGjGgftHPPPVcvvPDCsI/jeZ5yeaMg8GHA8zydetppuvTSS/X9731v4PFnn31Wn/3c5/TAAw/os5dfLklqa2vTRYsXa9OmTfrB97+v2tratF+3d/tGeW1NKexhSyqTwqWSPFnGk5qaZFpeTrsGAMhXbnvroH8nGwoyEQRSua5kSq77OOT0NkGm5ov+zW9+o9tuv13vrV+v6upqrV69WhdceKG+853v6Na/+Zsjtt+3b5++8MUvaskjj2js2LFpvWbXU4+oZ9Urw6wcAJCOoW4fxOPxjEw0lMp15YUXXtAv/vmfh3WbQEp9HZxMynkHwkz4whe+oJqaGt13333atWuXPnP55brqi188ahCQpLFjx+rEuXO1/KWX0n7NoAyFBIB81N9SsO++uwZ1NMzU7+ZUrytBl9MwkKk3LRQK6fbbbtN9v/iFPnP55Tr1lFN0zz33DNpm3759am9vlyS1t7dr+UsvadbMmRl5fQBAbhwRCjI0kiCZ60qm5fJDZk5vE3iep2iG1pfu7OzUpMmTNWPGDD2/bJnKygYP21u5cqVuuukmGSW+4Td85Sv66le/mv7r/fFh9a5+dZhVAwAyJTJxusoXXyG7rj4jxzvWdUWSLvvkJ/X222+rq6tL1dXV+vclSzRv3ry0Xi9SVDTk+jfZltMOhJm8N3LrbbdJkpqbmo46ecNpp52mFStWZOz1AAD+EJk4XZXnX6aiGXMUj8czNpHdsa4rkvSH3/8+I68lZfaamKq86DNw99136+mnn9bzy5Yp7rp68MEHs/6auV9jCgAKW2TidNX+1TdU95W/VXHDXFmWlbHfzbm4ruRSzicd6u3pGdb+DzzwgG7/5jf11FNPaf6ZZ+rHP/6xfn3//Xp3zRqFw9lb8a/35WfUt/GdrB0fAPJZa1md4qGIiqNdKo52K+RGB13I4y37ZfqOPtrs8JaAD36adl132KsV5uq6UlxSkrVjH0vOw0C0ry/t8ZXPPPOMrrjySj304IP69Kc/LSnROfC444/Xj374Q11zzTUZrHSwXA4BAYCg297arWVbmgf+7ViWyiKOSiOOSsOOrA1vydmzWcW97Srq6VBJd5tKx45X1RAhoN9wh6zn6rpi27YiRUVZOXYych4G4rFYWssXr1y5UhctXqy77rpLt9x886Dnvvvd7+rRxx7T22+9lbXFH3KZ4AAg6FzP6NHVe9TnJv9hMOJYqigKqaIorMqikCqKQ4m/i0IqDtkDASHdFudcXldCoZBCWWx1OJachwHPdRUdZpPOSMt1ggOAfPDazlat29+ZkWOFbGsgGMypLVFFJFhd4iKRiOwcrlyY8+9WsgtP+EkuF5MAgHzRUJO5lVvjnlFLT0yjSsIaXV6cseOOlFxfC3N+JbYsS6GAXVxzNQ4UAPKFZ4w6o3E5Gex7dfL4Sp1SXxW439GhUCjnfdB8cRV2QqG0+g3kgh/eNAAIqs6+uDY2dWljc6d6YplbnKc/CEjvf8gMynXFD63Nua9AiTfOcRy5GZpGMpv88KYBQJB4xmjXwR5tONCl3e3DX5zugw4PAv2C8iHTcRxffMD0zZUtFA77PgyEwmFfvGkAEATZagU43CnjK3XyB4KAdKh1IBzOyAqG2ZTLEQSH800YsCxL4Uhk2JNFZItl21kbpggA+SLbrQCHO6W+UiePPzII9OtvcTZpzmWTbeFIxDcfMH0TBqRDb5xtpz0JUTZFaBUAgCHFXE+bmrv07r4OdUWz38p7rCAgJT5kRsJh9WVoQbxMsn32ATPn8wx8kDFG0b6+nC7l+EHhcJi+AgBwFL1xV+v3d2r9/s6UJhAajqFuDQzFjccV89HtAsuyFCkq8tUHTN9d4SzLUiQS8U2SCxEEAOAInX1xvbuvQxubuuSO4Ie32WMrdNL4ypT2cUIhGck3/QciPro90M93LQP9PM9TLBrNaQtBKBwO3BwIAJBNLd1RrdnXoW0t3Rrp384NNWVaMKU67QtpPB7PaSDo7xvnx3kQfBsGpNzeMghHIr66nwMAuWKMUWNnn9Y0dmhPljsFDmXyqBKdO71G9jA/UWdiVcN0+PHWwOF8HQakxA+hG4+P2HhRy7YVDod9mdwAYCR5xmhHW4/WNHaouTt7F9CQbWlqdamm15TquU1NinuDL0vjKop0QUOdHDszF1LP8xSLxUZslEEoFJLj8wnrfN8G3j9W1HacrL95oXDYNxNAAECuuJ7R5kMjA9r7svdBrLokrFm1ZZpeU6aIk/gANqW6RJubuwe2qSmN6KMzajMWBKRDi81FInJdN6u3DYL04dL3YaBf/5vneZ7i8XhGQ0EQUhsAZFs07um9A51au79DvfHsfPDqbwWYVVem2tIjO9I11JQNhIGq4pAumFmrsJP5i2n/lMWO42S89dmybYVCIdm2HZjrSmDCgPT+tMWO48jzPLmuKzfNN9C2bTkBe7MAIBu6onGt29+pDQc6FfOyc+f4aK0ARzO2vEgVRSF5xujCmXUqDmW371Z/67MTCiWuK/F42nPdOIfCRRBaAj7I930GkmGMked5Mp6X+PuD/yXLkm1ZiQv/oYs/AQBAoWvriendfR3a0tKlbGQAS9K00aU6fkz5UVsBhrLhQKfGVhSpqjg3U/UaYxJ/Dl1TPGOkD1xXrMOuKfnwoTIvwgAAIHkdfXGt3N2mba09WTl+yLY0s7ZMs8dUqLwoUA3QBYt3CQAKRF/c0zuN7Vq3vyMrLQFFjq0TxpTruDHlWW/eR2YRBgAgz7me0XsHOrV6b3tWpgwujziaPbZCM2vLFArg/XIQBgAgb5lD8wS8ufugOrIwRLC6JKy54yo0tbp02JMBIbcIAwCQhw509emNXW3a35n5yYLGVRRp7tgK1VcWB77jHBIIAwCQRzr74npz90Fta+0+9sYpmjKqRHPHVai2rCjjx0ZuEQYAIA9E455WZ6FzoG0lJgKaPbYiZ0P9kH2EAQAIMM8kOgeu2pPZzoFhx9JxdeWaPaZCJWFGBuQ7wgAABFB/58CVuw9mdP2AkrCt2WMqNKuu/ENnCkR+IQwAQMA0dfXp9V0Htb+zL2PHLAnZOqW+SjNqyjK6KBCCgTAAAAHR2RfXyj0HtbUlc50DQ7alOWMrNGdsRVYWBEIwEAYAwOeirqd39rZrbYY7BzbUlOnU+iqVRugTUOgIAwDgY7sP9uiV7a3qirkZO+b4iiLNmzhKo0sjGTsmgo0wAAA+1Bf39PquVm1uztwtgVHFIc2bOIrJgnAEwgAA+MyOth69uqNFPbHMDBUsCdk6ZUKVGmrKmDYYR0UYAACf6I27em1nW8Y6CDqWpbnj6ByIYyMMAIAPbGvt1oodreqNZ6Y1oKGmTKfUV6oswq95HBs/JQCQQz0xVyt2tGp7W09GjkfnQKSDMAAAOWCM0ZaWbr2+sy0j0whXHeocOIHOgUgDYQAARlhXNK5Xd7Rq18HeYR+rOGTr1PoqNdTSORDpIwwAwAgxxmhTc5de39WmmDu82YNsS5oztlInjqNzIIaPMAAAI6CzL65XdrRqT/vwWwNqSsNaOHW0qkvoF4DMIAwAQBYZY/ReU5fe3NWm+DDnErYt6dT6Ks0eW8EtAWQUYQAA0tDV1aX9+/errKxMY8aMOeo27X1xvbKtRY0ZWF2wriyihVNHq6o4POxjAR/EjSYASFJjY6NuvPFGzZkzRxUVFZo+fbrGjh2ru+66a9B2xhit3deh37/bOOwg4FiWzpg4ShcfN4YggKyxjDEZXAMLAPLXTTfdpMcff1yf+9zndPrpp2vq1Kl66KGH9OKLL2rr1q2SErMILt/aot0Z6BswtrxIC6aOVmURjbjILloGACBJe/fu1bx583TLLbcoEokMBALXTawo2NQV1X+t2zfsIBCyLc2fXK3Fs+oIAhgR/JQBQJJOPPFEfe9739NTTz0lSdq0adPAc+8d6NRrO1s1zD6Cqq8s1tmTq1VOCMAI4qcNAJJ055136owzztD69ev17W9/e+DxvrinV3e0DuvYYSfRN6ChpowZBDHiCAMAkKRQKKRPfvKTqqioGPT4cIcMTqwq1lmTq1lUCDnDTx4A5EiRY+vMSaM0bXQprQHIKcIAAKTAM0br93cM+ziTR5XorMnVKgk7GagKGB7CAAAkqTvq6uk127Vy/aZBj7vxmPbt3KYxE6cc8xN+ccjW/MnVmlpdms1SgZQQBgAgCY0dvfrXp1/Undddrt6uzkHPHWw+oG9/5iM67rSz9Lf/378PGQgmVhVr4dTRKg7RGgB/IQwAwIcwxmjNvg69tfugVvzlWZUVF+nZpxNDCydMmKCvfOUruuiii/TGG2/o1ltvVev+Ro0eO37QMSxJp9RX6cRxFfQNgC8RBgBgCNG4p+XbmrXzYGISobLKKvX09OiBBx6QpIG/JWnPnj1ynJCKSkoGHaM4ZOsj02o0vrJ45AoHUkQYAICjaO6OatnmJnVG3YHHzrr4U9q+fo2ef/3tI7a3bVt/9e2/V1nlqIHH6soiOnd6DUMG4XusTQAAH7CxqVOv7hjebIInjCnX6RNGybG5LQD/I64CwCFxz9OKHW3a1NyV9jFCtqWFU0Zr6mhGCyA4CAMAIKkrGtdzm5rU0hNL+xhVxSGdP6OWpYYROIQBAAWvrSemP286oK7D+gekalp1qc6eUq2ww2KwCB7CAICCtr+zT0s3NSnqemntb1vSGRNH6bi6coYNIrAIAwAK1o62Hr2wpVlumv2oy8KOzp1Ro7qyogxXBowswgCAgrThQGLEQLoDBuori3XONGYTRH4gDAAoKMYYrdrbrlV729M+xsnjK3XS+ErZ3BZAniAMACgYnjFasaNVG5rSGzpY5Ng6Z9poTagqOfbGQIAQBgAUhLjn6YUt708tnKqa0ojOm16j8iJ+bSL/8FMNIO/1xV0t3dSkA13RtPafUVOqsyePZjZB5C3CAIC81hmN688bD+hgbzyt/U8cV6lT6ysZNoi8RhgAkLdae6J6dmOTemLpTSY0f9IoHT+mIsNVAf5DGACQlxo7evXc5ibF3NQHD9qWdM60Gk2tZn0BFAbCAIC8s621Wy9ubU5r1cGwY+mjM2o1rqI484UBPkUYAJBX1u3v0Gs729LatzTs6IKZtaouiWS2KMDnCAMA8oIxRm/tOah3GjvS2r+qOKQLZtapPMKvRRQefuoBBJ5njF7e3qLNzd1p7V9XFtHHGmpVxNTCKFCEAQCB5npGy7Y0aVeakwlNqirRR6aPVshm6WEULsIAgMDyjNELW5vTDgKzass0f3I1awyg4BEGAASSZ4yWb23RjraetPY/5dBiQ0wmBBAGAASQMUavbG/V1tbU+whYks6aXK1ZdeWZLwwIKMIAgEAxxmjFzjZtak595UHHsvSR6TWaPIpVB4HDEQYABIYxRm/sPqj3DnSmvG+RY+ujDbUaU16UhcqAYCMMAAiMVXvbtXZf6vMIlEUcXdBQp1El4SxUBQQfYQBAILzT2K5Ve9tT3m9UcVgXzKxVGZMJAUPi7ADge2v3dWjl7oMp71dVHNLiWXUqDjOZEPBhmGUDgK9tONCp13e1pbxfRVFIF80aQxAAkkAYAOBbm5u79MqO1pT3K4s4WjyrTqUEASAphAEAvrSttVsvbWtJeb+SsKPFs8bQRwBIAWEAgO/sbOvRC1uaZVLcrzhka/GsOlUUEQSAVBAGAPjK7vZeLdvSlHIQKHJsXTSrTlXFDB8EUkUYAOAbjR29+sumJnkpJoGwbenCmXWqLolkpzAgzxEGAPjC/s4+Ld3UJNeklgRCtqULZtappowgAKSLMAAg55q7o/rzpgOKp9gk4FiWPsYUw8CwEQYA5FRrT0zPbjigmJtaELAt6fwZNRpXUZylyoDCQRgAkDM9MVdLNx1Qn+ultJ8l6dzptZpQxeqDQCYQBgDkhOsZLdvcpK6om9J+lqRzprEMMZBJhAEAI84Yo1e2t2h/VzTlfRdMHa1po0uzUBVQuAgDAEbcu/s6tLmlO+X9zppcrYaasixUBBQ2wgCAEbWjrUdvprEC4RkTR+m4uvIsVASAMABgxLR0R/Xi1uaU9zu1vkqzx1ZkoSIAEmEAwAjpibl6bnNTynMJnDiuQieNr8xSVQAkwgCAEZDuyIGp1aU6tb4qS1UB6EcYAJBV6Y4cqCmNaOHUalmWlaXKAPQjDADIqnRGDpSGHX20oVYhm19RwEjgTAOQNemMHHAsSx9tqFVp2MlSVQA+iDAAICvSHTlwzrTRqillBUJgJBEGAGRcuiMHTqmv1JRqZhcERhphAEBGuZ7RX9IcOXDSOIYQArlAGACQMf0jBw4wcgAIFMIAgIxh5AAQTJx9ADKCkQNAcBEGAAwbIweAYCMMABiWXkYOAIFHGACQNmOMlm9rSXnkwDRGDgC+QhgAkLY1+zq0u703pX1qSiNawMgBwFcIAwDSsq+zT2+l2GGQkQOAP3FGAkhZb9zVC1ualUovAUYOAP5FGACQEmOMlm9tUXcstX4CjBwA/IswACAl6fQTYOQA4G+EAQBJ259GP4GJVcWMHAB8jjAAICm9cVfPp9hPoCziaNHU0YwcAHwulOsCAIw848blHTwgE+2V3HjiQSckK1Isu6pOljP4V0M6/QQsSedOq1FRiA6DgN8RBoAC4HV3KL5znbyWRrnNu2XamyQzxGd8y5JVWSunZoLs0eMUmnSC9sZCKfcTOH1ilerKizJQPYBss4wZ6jcCgCAznid372bFNr4pd/fGxIOWJRkvuQNY9kBgcCbO1J6JZ2hFR5Fi7rF/ZUysKtZHZ9RyewAICMIAkGeMMYpvXa3o28/J9HQcuqgnGQCGcugY3eVj9caUC3XADQ+5aVnE0WUnjOX2ABAghAEgj3idbepb8Qe5jVuz9xqytGHKuVpbMu2IzoSWpIuPG6Mx3B4AAoUwAOQBY4xi772m6NtLJc8bfkvAsVi2mkrH6vXJF6jbvN/16PQJVZrLMEIgcAgDQMAZz1Pfa/+l+Oa3R/y1o3ZEbzdcop12Nf0EgAAjDAABZjxXvS8+JnfXe7mrQdKuqYs09YxzVBIZui8BAP8iDAABZTxPvS89IXfH2lyXIklyJs9W8cLLZbEiIRA4nLVAQMXWvuSbICBJ7o61iq19KddlAEgDYQAIILdlr6Krl+W6jCNEVy+T27I312UASBFhAAgY48bVu/yJXJcxpN7lT8j0T3EMIBAIA0DARFcvk+loGXo64VwyRqajxZetFgCGRhgAAsREexVbv0JKae3AkZaY88BEU1vLAEDuEAaAAIltflvykl85MGfcuGJb3s51FQCSRBgAAsJ4nmLrX811GUmLrXtVxsvyTIgAMoIwAASEu3ezTHf7iL1ed19Us7/xD/ofDz+T1v6mu13u3s0ZrgpANhAGgIBwG7ckVg8cIT/9/1/QvBkT0z+AZWd1wSQAmUMYAALCbdqd/QWIDtnU2KwNe5t00Skz0z+I8eQ27cpcUQCyhjAABIDxPHnDnMzH8zyd9t//r/7nkj8NevzPqzdp9DXf1e9WvDvw2N89/IzuuuKCYb2eJHmtjfQbAAKAMAAEgHfwwLBHEdi2rdsvO0e//vMbau3qkSS9s71RX/7Zf+jOKz6mz8yfI0n6rzfXq2FcjWaOrx123XLj8tqbhn8cAFkVOvYmAHLNdLVl5DhXLjxRP/rdMv3zM6/qy+eeps/979/qC4tO0jc+sXBgm9c37dTjr67Rf762Vp29UcVdV5UlRfrbz5yXXu2drdKoMRmpH0B2sGohEACx7e+qb/njGTnW/Utf112PLlV9daUm147SI7d+Qc4QKw3+5oW3tHbXfv3gqsVpv17Ros8qPGVO2vsDyD5uEwBBkMH77lcsPEk9fTEZY/TAzZ8bMghkDH0GAN/jNgEQAJaTuVP19of+KElq7uiWY1sfuu3VHzl12K+XydoBZActA0AQhCMZOczfP7pUz7y9Qc/ddYPinqd/fX5lRo77oTJUO4DsIQwAAWCPGjvsYzz4lzf1f596Rf9x+1U6cco43XzxWbr3v15SLJ7dtQ4yUTuA7CIMAAFgl5TLKi5Le/8/rdqo2x96Ur+86XKd2TBJkvTfLpyvjp4+PfLSqkyVeQSruFx2SXnWjg8gMwgDQEDYtRMlffg9/qN5a+seffln/6G//8KF+tQZswcerywt1n+7aL7+4Q/L5Wajk59lHaoZgN8xtBAIiOia5Yqu/osUlFPWshQ5+XxF5izKdSUAjoGWASAgQlNmBycISJIxCk2efeztAOQcYQAICLtitJz6BslK/VbBiLNsOfUNsitG57oSAEkgDAABEj7uzGC0DhhP4ePn57oKAEkiDAAB4oyfIau8Wul0JBw5lqzyajnjpue6EABJIgwAAWJZlormXSzJz60DRkVnfFxWEG5nAJBEGAACJzRhpkINp/m074ClUMPpCtU35LoQACkgDAABVHTaRbJKK/0VCCxLVlmVik6/MNeVAEgRYQAIICscUfHCy+WvvgOWihddLivEWgRA0BAGgIBy6iap+COf90frgGWp+COfl8OMg0AgEQaAAAtNPE7F53xesmzlppXAkixbxed8XqGJx+Xg9QFkAtMRA3nA3bdNPcuWSG5s5OYhsCzJCavkvC/KGTtlZF4TQFYQBoA84XV3qO+1J+Xu3jAir+dMPE5FZ1wiu7RiRF4PQPYQBoA8YoyRu3Odel97Uor2Zr6VwLKkSLGKz7xUocknZPbYAHKGMADkIdPXo9iG1xXb8IZMb2fiIp7uqX5oX6u4XOFZ8xSedYasopLMFgwgpwgDQB4znid37ybFNr4pd/cmDcxcaNmS8Y6+06DnLDkTGhSeOS8xFbJNn2MgHxEGgAJh4lF5rfvktTbKbWmU17xHJtoj48YlSZYTkhUpkV1TL2f0ONnV42RXj2XeAKAAEAYAAChwtPkBAFDgCAMAABQ4wgAAAAWOMAAAQIEjDAAAUOAIAwAAFDjCAAAABY4wAABAgSMMAABQ4AgDAAAUOMIAAAAFjjAAAECBIwwAAFDgCAMAABQ4wgAAAAUulOsCDhc59XrZoYgs25FlO3LC739t2fb7zzmO7FBE9sBzzhHPWbYj27Zk2ZYcx5b1ga9t25LtWAPbfOhzliUnZMuxLTm2pcihr0MD/3bef855f7vQYds6R/vasmRblhxLCjv2wNchx5ZjKfFv21LYto7ydeL5sG0PfO1YlixLsi3JsnTo+JIlybEt2VLi/2Jr4Gvbkhzr8K8Tx7CMkYwny4tLg772En+8oZ+zjCe57vtfe3HJc2U8T4pHZVxX8rzEY/GYjOcmvo7FpP6v+7ft3y4WfX8fz5UXi8u4noznyYvG5bmJfYzryYvF5bnvf20Ofe3G4jKHbedG44d97cp4Rp5rDv370P6eSTznGhnXyHM9uTHv0DGN3Jh7aJ/39/OMkWuMop6Ra/SBrz/478TXnhJfu0aHnnv/638223J6XmYK5zfnN+e3f89vWgYAAChwhAEAAAocYQAAgAJHGAAAoMARBgAAKHCEAQAAChxhAACAAkcYAACgwBEGAAAocIQBAAAKHGEAAIACRxgAAKDAEQYAAChwhAEAAAocYQAAgAJHGAAAoMARBgAAKHCEAQAAChxhAACAAkcYAACgwBEGAAAocIQBAAAKHGEAAIACRxgAAKDAEQYAACh0Jk/19vaaO++80/T29ua6lCP4uTZjqG84/FxbPvHz99nPtRlDfcPh59qGyzLGmFwHkmxob29XVVWVDh48qMrKylyXM4ifa5Oobzj8XFs+8fP32c+1SdQ3HH6ubbi4TQAAQIEjDAAAUOAIAwAAFLi8DQNFRUW68847VVRUlOtSjuDn2iTqGw4/15ZP/Px99nNtEvUNh59rG6687UAIAACSk7ctAwAAIDmEAQAAChxhAACAApd3YeCb3/ymzjnnHH3pS19SNBod9FxPT48uvfRSnXvuubrwwgvV0tLiq/r6/fCHP9S8efNyXlM8Hte1116rc845R9/4xjdGrJ5k6+s30t+vww1Vmx9+1vIR53fmauL8PrZCOr/zKgy89dZbamxs1IsvvqjZs2frscceG/T8U089pblz5+r555/XFVdcoX/7t3/zVX2S1NHRoTVr1viipj/84Q+aOHGiXnzxRXV3d+vll18esbqSqU8a+e9XsrXl+mctH3F+Z7Ymzu/0a8v1z1o25FUYeOWVV3TRRRdJki6++OIjfrhnzpyp7u5uSVJbW5vq6up8VZ8k/eM//qNuvvlmX9SUTL25rE8a+e/X4T6stlz/rOUjzu/M1sT5/eEK7fwO5bqATGpra1N9fb0kqaqq6oimmxkzZmjNmjWaO3euLMvSihUrfFXfwYMH9c477+iOO+7wRU1tbW0D828frd5c15eL71eyteX6Zy0fcX5ntibO7/Rry/XPWjYEsmWgsbFRixYtOuKPMUbt7e2SEm/k6NGjB+330EMP6bzzztOaNWt0991367vf/a6v6rv33nt1yy23ZKWmoVRXVw9Z04c954f6cvH9OtyH1TZSP2v5iPM7czi/01do53cgw8C4ceO0fPnyI/5ccskl+tOf/iRJeuaZZ7Rw4cIj9u1/Q0eNGqW2tjZf1bdp0yZ9//vf18UXX6yNGzfqRz/6UVbqO9xZZ501ZE0f9txI+bAacvH9SrY2aWR+1vIR53fmcH5npzYpD8/v3K2enB233367WbRokbnqqqtMX1+fMcaYr371q8YYYw4ePGguueQSc+6555qFCxea9957z1f1He7000/PWU399cRiMfPlL3/ZLFq0yHz9618fsXqSre9wI/n9OtxQtfnhZy0fcX4PvybO7+QV0vnNdMQAABS4QN4mAAAAmUMYAACgwBEGAAAocIQBAAAKHGGgADz44IMaNWpURo61bds2WZalUCik3bt3D3pu7969CoVCsixL27ZtG/Tc448/rvPOO09VVVUqLy/XSSedpO9+97sDE3lkskag0Fx77bWyLEs33njjEc997Wtfk2VZuvbaawcea2xs1Ne//nVNnz5dRUVFmjRpki677DItXbp0YJupU6fq3nvvHYHq4QeEAaSlvr5e//qv/zrosYceekgTJkw4Ytu/+7u/05VXXqkzzjhDTz31lNasWaN77rlHq1atyos5vQE/mDRpkpYsWaKenp6Bx3p7e/XII49o8uTJA49t27ZNp59+up577jn95Cc/0TvvvKOnn35a559/fs6m/kXuEQYC4Omnn9aiRYs0atQo1dTU6NJLL9XmzZslScuWLZNlWYMmvXj77bcHPp0vW7ZM1113nQ4ePCjLsmRZlu666y5JUmtrq7785S+rurpapaWl+vjHP66NGzcmVdM111yjBx54YNBjDz74oK655ppBj7322mv6wQ9+oHvuuUc//elPtWDBAk2dOlUXXnihHn/88SO2B5Ce0047TZMnT9YTTzwx8NgTTzyhSZMm6dRTTx14rL+l4LXXXtPnPvc5zZo1S3PmzNFtt92mV199NRelwwcIAwHQ1dWl2267Ta+//rqWLl0q27b1mc98Rp7nHXPfBQsW6N5771VlZaX27t2rvXv36pvf/KakRNPiG2+8od///vd65ZVXZIzRJZdcolgsdszjfvKTn1Rra6uWL18uSVq+fLlaWlp02WWXDdrut7/9rcrLy/W1r33tqMfh1gCQOdddd92gkH7//ffr+uuvH/h3S0uLnn76ad18880qKys7Yn/Ox8KVVwsV5avPfvazg/7961//WmPGjNHatWuPuW8kElFVVZUsy9K4ceMGHt+4caN+//vf66WXXtKCBQskJS7ckyZN0n/+53/q85///IceNxwO6+qrr9b999+vRYsW6f7779fVV1+tcDg8aLuNGzdq+vTpRzwOIPP+6q/+St/5zncG+va89NJLWrJkiZYtWyYpMcWvMUbHH398bguF79AyEACbN2/WVVddpenTp6uyslLTpk2TJO3YsSPtY65bt06hUEjz588feKympkbHHXec1q1bJ0n6+Mc/rvLycpWXl2vOnDlHHOOv//qv9eijj6qxsVGPPvrooE8g/Ywxsiwr7ToBJK+2tlaf+MQn9NBDD+mBBx7QJz7xCdXW1g483z/hLOckPoiWgQC47LLLNGnSJP3yl79UfX29PM/T3LlzFY1GVV5eLun9k1xSUs38Q81CffjF+1e/+tVAZ6SjfbKfO3eujj/+eH3xi1/UCSecoLlz5+rtt98etM2sWbO0fPlyxWIxWgeAEXD99dcPrPb385//fNBzM2fOlGVZWrdunT796U/noDr4FS0DPtfc3Kx169bpjjvu0Mc+9jGdcMIJam1tHXi+rq5OUmJYX78PXpAjkYhc1x302OzZsxWPxwetw93c3KwNGzbohBNOkCRNmDBBDQ0Namho0JQpU45a3/XXX69ly5YdtVVAkq666ip1dnbqvvvuO+rzebHaF+AjF198saLRqKLRqBYvXjzoudGjR2vx4sX6+c9/rq6uriP25XwsXIQBn6uurlZNTY3+5V/+RZs2bdJzzz2n2267beD5hoYGTZo0SXfddZc2bNigJ598Uvfcc8+gY0ydOlWdnZ1aunSpmpqa1N3drZkzZ+pTn/qUbrjhBi1fvlyrVq3S1VdfrQkTJuhTn/pU0vXdcMMNOnDggL7yla8c9fn58+frW9/6lm6//XZ961vf0iuvvKLt27dr6dKl+vznP6+HHnoovW8MgKNyHEfr1q3TunXr5DjOEc/fd999cl1XZ555ph5//HFt3LhR69at089+9jOdffbZOagYfkAY8DnbtrVkyRK9+eabmjt3rm699Vb99Kc/HXg+HA7rkUce0fr163XyySfrxz/+sb73ve8NOsaCBQt044036sorr1RdXZ1+8pOfSJIeeOABnX766br00kt19tlnyxijP/7xjyk154dCIdXW1ioUGvqO049//GM9/PDDWrFihRYvXjwwjOmkk05iaCGQBZWVlaqsrDzqc9OmTdPKlSt1/vnn6/bbb9fcuXN14YUXaunSpfrFL34xwpXCL1jCGACAAkfLAAAABY4wAABAgSMMAABQ4AgDAAAUOMIAAAAFjjAAAECBIwwAAFDgCAMAABQ4wgAAAAWOMAAAQIEjDAAAUOD+H09tb3u7E2wvAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pc_alpha_bootstrap = 0.1\n",
    "boot_samples = 200\n",
    "\n",
    "# The block-length of the bootstrap can optionally be used to better deal with autocorrelation, \n",
    "# but its effect was not yet evaluated.\n",
    "boot_blocklength = 1\n",
    "\n",
    "## Create PCMCI object to call run_bootstrap_of\n",
    "pcmci = PCMCI(dataframe=dataframe,\n",
    "        cond_ind_test=ParCorr(),\n",
    "        verbosity=0,\n",
    "        )\n",
    "\n",
    "# Call bootstrap for the chosen method (here 'run_pcmciplus') and pass method arguments  \n",
    "results = pcmci.run_bootstrap_of(\n",
    "        method='run_pcmciplus', \n",
    "        method_args={'tau_max':tau_max, 'pc_alpha':pc_alpha_bootstrap}, \n",
    "        boot_samples=boot_samples,\n",
    "        boot_blocklength=boot_blocklength,\n",
    "        seed=123)\n",
    "\n",
    "# Output graph, link frequencies (confidence measure), and mean test statistic values (val_mat)\n",
    "boot_linkfreq = results['summary_results']['link_frequency']\n",
    "boot_graph = results['summary_results']['most_frequent_links']\n",
    "val_mat = results['summary_results']['val_matrix_mean']\n",
    "\n",
    "# Plot\n",
    "tp.plot_graph(\n",
    "    graph = boot_graph,\n",
    "    val_matrix= val_mat,\n",
    "    link_width = boot_linkfreq,\n",
    "    var_names=dataframe.var_names,\n",
    "    ); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we can see that Bagged-PCMCI+ gets the correct graph (link-wise majority vote among the $B$ resamples) and also provides a useful estimate of the confidence for links via the widths of the arrows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap aggregation with optimized ``pc_alpha``\n",
    "\n",
    "For a range of algorithms and conditional independence tests, it is also possible to set ``pc_alpha`` to a specficied list or to None (then the list [0.001, 0.005, 0.01, 0.025, 0.05] will be used). ``pc_alpha`` will then be optimized across the values specified in the list using the score computed in cond_ind_test.get_model_selection_criterion(). This is not possible for all conditional independence tests. \n",
    "\n",
    "This approach can also be combined with bootstrap-bagging. ``pc_alpha`` is internally and individually optimized for each of the  $B$ bootstrap resamples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAHBCAYAAAD0E7h1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA94klEQVR4nO3deZRcdZ3//9e9t6p63ztJJ53OnpCNQEjYkiCIIqCgiAuKaJBRvyr6GwXGc/zqjMDXFYcRmRFnRmWZUcBhcY6jgjJgAgQIS1gSaLKSvTtJb+m9q+rez++PTjdpslVVV3XdW/V8nJOTTtW9t97p6tv3VZ/7WSxjjBEAAMhbdrYLAAAA2UUYAAAgzxEGAADIc4QBAADyHGEAAIA8RxgAACDPEQYAAMhzhAEAAPIcYQAAgDxHGAAAIM8RBgAAyHOEAQAA8hxhAACAPEcYAAAgzxEGAADIc4QBAADyHGEAAIA8RxgAACDPEQYAAMhzhAEAAPIcYQAAgDxHGAAAIM8RBgAAyHOEAQAA8lwo2wUAGFvGmME/nicz+MDgE5YlS5Jl27IsS5ZlZbFKAGOJMADkOGOMPM+T57ryPE9m6OJ/ApZlybZt2Y4j+1BAAJCbLJPobwYAgWKMkRuPy3XdhAPAsViWJcdx5IRChAIgBxEGgBxjjFE8Hpcbj2fk+E4opBChAMgphAEgh7iuq1g0OiavFY5E5DjOmLwWgMwiDAA5wBijWCwmz3XH9HVtx1E4HKaVAAg4wgAQcMYYRQcGRt0vIFWWZSlSUEAgAAKMeQaAAMt2EPBLDQBGhzAABJSfLsJ+qgVA8ggDQEDFYjFfXXyH+i0ACB7CABBAruuOeWfBRHiuK9eHdQE4PsIAEDDGmDEbPpiKWDTqqxYLACdGGAACxs9BYEgQagTwNsIAECBD6wz4XTJrIADIPsIAECDxDE0xnAlBqhXId4QBICCGFh4KCjcep3UACAjCABAQYzV6YNfu3XrfhRdq8Wmn6fQzztBDDz+c8rH8OOIBwJGYjhgIiFg0OibD9pqamrR//36dcsop2r9/v85etkyvvfqqSkpKkj6W4zgKRyIZqBJAOoWyXQCAxIxVx8GJEydq4sSJkqTx48erurpabe3tKYWBIHR2BMBtAiAQjDGjvv/ueZ5OOfVUfevb3x7x+GOPPabyioqj3g546aWX5HmeGiZPTuk101E3gMwjDAABkI4Lqm3b+rsbbtAvfvELtbe3S5Jee+01feqqq3TTTTfpI5dfPmL71tZWfe7zn9e//Mu/jOp1CQOA/9FnAAgA13XTMpFPPB7XyYsW6apPfUorV67Uueedp0svuUS33XbbiO0GBgb0gUsu0TWf/ayuvPLKUb1mOBKR4zijOgaAzCIMAAHgxuNpWwTol7/8pb5z442aNGmSpk6Zot/+9rcjLtbGGK28+mrNmT1b337HLYVUhMNhOSG6JwF+RhgAAiCdYaC7u1sNU6Zo5syZWr1q1REdA9c884wuuOACnbxw4fBjv/rVr7TwsH8ngzAA+B9nKBAElpW2Q339uuskSa0tLUdtvl++bJl6e3rS9nrprB1AZtCBEMgjN910kx599FGtXrVKcdfV3Xffne2SAPgAYQAIANse/al611136ae3364HH3xQixYt0leuvVb/9JOfpO32w7Gko3YAmcVZCgSAZVmjam7/85//rK99/eu681e/0plnnCFJ+tKXvqSuri7de++96SrzSJY1WDsAXyMMAAGR6ifsdevW6VNXXaXvfe97uuyyy4YfLy8v15e++EX94623ZmyaY1oFgGBgNAEQEPF4XPEMN+mnWygcVoiRBIDvEduBgAjixD1BrBnIR4QBICAsywpUs7tt2/QXAAIiOL9ZACgUDme7hIQFqVYg3xEGgAAJyqftoLViAPmOsxUImHAkku0STigINQJ4G2EACBjbtn3dQz8UCtEqAAQMZywQQE4o5MvbBZZlsSgREECEASCALMvyZVN8OBLxZUgBcHyEASCgbNtWpKAg22UMixQUcHsACCjOXCDAbNtWqP+g1N2a1ToIAkCwcXMPCDAz0Ctvywuy41GZulkyE2ZJ1thdlIduVxAEgGBjbQIgoIzxFH/9SZmuw1oFSqrkTTlFKijO+Os7oZBCPu3ICCA5hAEgoOI7N8jbs/HIJ5ywTMNCmcqJGXldWgOA3EMYAALIa29W/M01x93Gqp0ir36+TJpuG1iH5jcIyiyIABJHnwEgYMxAn+JbXjjxdi075YQjsqecLNd15bqulGz2tyw5jiPHcWgJAHIYYQAIEGM8xTc/L8WjJ964oETO5PmybFu2bSscDssYI+N58g79bYzRUDywNHgLwLJt2Yf+pgUAyA+EASBA3F1vyHS1nHhDy1ZozpmyQiNXDrQsS5bjMKYYwAj8TgACwutoPnqHwaNwpp4su7QqwxUByBWEASAAzECf4ptP3E9Akqzqetl1MzNcEYBcQhgAfM54ruKbnku4n0Bo5hLu9QNICmEA8DFjjNyt62S620688TH6CQDAiRAGAB/z9m6S17IzoW3pJwAgVYQBwKe8tj1yd25IaFurehL9BACkjDAA+JDX05Fwh0EVFNNPAMCoEAYAnzHRfsXffEby3BNvbNkKzTlLViiS+cIA5CzCAOAjxnMV3/isFO1LaHtn1lL6CQAYNcIA4BODIwdeSmzkgCR78jw5tQ0ZrgpAPiAMAD7h7dkor2VXQtvaNZPlTJ6X4YoA5AvCAOADXuseubteT2hbq6RKDh0GAaQRYQDIMq+7PaEliSVJ4UKF5p4ty2GNMQDpQxgAsshE+wY7DCYycsB2FJq7TFakKPOFAcgrhAEgS4zrKv5m4iMHQowcAJAhhAEgCwZHDrwo09Oe0PZOw3zZNZMzXBWAfEUYALLA2/2mvNbdCW1r10yWXT83wxUByGeEAWCMuS275e5+I6FtrdIqOTOXMnIAQEYRBoAx5HW3yU105ECkSKGTzpblOJktCkDeIwwAY8QM9A12GDTeiTdm5ACAMUQYAMaAceOKb3xGivUntH1o1umySyozWxQAHEIYADLMuK7iG5+R6elIaHtnygLZNfWZLQoADkMYADLIeJ7im56TOXggoe3t2imyJ52U4aoAYCTCAJAhxvMU37xWpqM5oe2t0mo5M09j5ACAMUcYADLAGCN3ywsybXsT2yFSPDhywGbkAICxRxgA0mxodsFEJxUaHDlwtqxIYWYLA4BjIAwAaWSMkbvtZXkHdia4h6XQ7DMYOQAgqwgDQJoYY+Ruf1Xe/rcS3seZfbrs6kkZrAoATowwAKSBMUbuzg3ymrcmvI8zc4mc2oYMVgUAiSEMAGng7m6Ut3dTwts70xfLGT8tcwUBQBIIA8AouXvelLe7MeHtnWmL5NTNyGBFAJAcwgAwCm7TZrk7X094e2fKQjkTZ2ewIgBIHmEASJHbvFXu9tcS3t6ePE9OPbMLAvAfwgCQAnf/drlvvZLw9vakOXImz8tcQQAwCoQBIEnugZ1yt76U8PZ23Uw5UxYyzTAA3yIMAEnwWvfI3fJiwtvbE6bLmXYKQQCArxEGgAR5bXsV37xWkkloe3vcVDnTFxMEAPgeYQBIgNexT/FNayWTYBComSxn5hKCAIBACGW7AMDvvNY9im9+XjJeQttb1ZPkzDqdIAAgMAgDwHG4TVvkbn814e2tyjqFZp8py6bRDUBwEAaAoxheayCJKYativEKnXQWQQBA4BAGgHcwnid364vyWnYlvI9VVqvQSWfLsp0MVgYAmUEYAA5j4jHFNz0rc/BAwvtYpdUKzVsmy+F0AhBM/PYCDjHRPsUb18j0Hkx4H6ukUqF5y2U54QxWBgCZRRgAJJm+TsUa10gDvQnvY5VWKzR3uaxQJIOVAUDmEQaQ97zOFsU3PivFownvY1VNUmj2GbIc+ggACD7CAPJasnMISJI9YYac6acyjwCAnEEYQN5ym7cmtfKgJDkNC2TXn0QQAJBTCAPIO8YYubtel7dnY+I7WZacGUvkjJ+aucIAIEsIA8grxvPkbntJ3oGdie9kOwqddJbsyrrMFQYAWUQYQN4wbkzxjc/JHNyf+E7hAoXmLpddWpW5wgAgywgDyAsm2q/4m2tkejoS36mwVOF5K2QVlmSsLgDwA8IAcp7p61Ks8ekU5hBYJitckMHKAMAfCAPIaanNITDx0BwCnB4A8gO/7ZCTjDHymrfI3bFeMibh/ewJ0w/NIcDKgwDyB2EAOcfEY4pvfUmmbU9S+zkN82XXz2UOAQB5hzCAnOL1dCi+6TmpvyeJvSw5M0+TM35apsoCAF8jDCAnGGPk7d8+OKNgElMLy3YUmnOW7CrmEACQvwgDCDzjxuW+9XJyEwlJzCEAAIcQBhBopq9T8Y1rZfo6k9uxsOTQHAKlmSkMAAKEMIDA2rn+BT33x4fU2tGhibXVWjJvpiZUV55wP6t8nEJzzmQOAQA4hDCAwNmw/jV9fuWn9dzLr0mSLMuSMUbhUEhf/ujF+n9fvFIFkfBR97Xr58ppmM+IAQA4DIOpESimv0f/cN1X1NLVq3vvvVfNzc3yPE+7du3SjTfdpH/778f09Z/ceeSOoYhCc5crNGUBQQAA3oGWAQSG196k+OYX1NHRoUWLFmlgYEArV67U9u3bdfrpp+tnP/uZSkpKdP311+lb13xU9eNqJB2aWnjOmbIKirP8PwAAf6JlAL5njKf4jvWKv/mM5Ma0YMYUPfzww7rmmmsUa2vS+SfP0K9//WvdfvvtuvLKK+W6nl54fbMkya6bpdCCcwkCAHActAzA10y0T/FNz8t0tQw/9oNrr9I5p87T3GkNmjutXnsOtOqePzwhx3HU3d0tSSooLFRozpmyayZnq3QACAzLmCQmbgfGkHdwv+Kbn5diA8fcpqunTxd+9UY1d/Vr/fr1uuWWW/TT236iPds2q2piwxhWCwDBRcsAfMcYI2/PRrm7Xj/udj19/brshu9rS1OLVq1apRdeeEE//vGP9Q9///cEAQBIAi0D8BXT1zW4yFBX63G36+0f0Idv+IFe2rxDf/nLXxSNRnXxxRfr/PPP1+9+9zuFQuRcAEgUvzHhC8bz5DVtlrvrjROuLdA3MKDLv/FDvbhpux599FEZY/SBD3xAfX19+tKXvqS2tjaNHz9+jCoHgOCjZQBZ5/V0yN3yokzvwYS2/+bP/lP/+rvH9Mgjj6ioqEgXXHCBDh58e99wOKwNGzZozpw5mSoZAHIKLQPIGuO5cnc1ytu7SVLimfSvL67XVVddpXe9613at2+f1q1bN+L5OXPm6IknniAMAECCCAPICq+zRfGtL0n93UnvO3F8rR555BFdfPHFRz+256m+vn60JQJA3uA2AcaUcWNyd2yQt29bSvtblXXa6pbrn27/Z7W1tR11myVLlujv/u7v6EQIAAkiDGDMeO3Nim9bJ0X7UtrfmbJA9qSTWFsAANKMj07IOBMbkLv9NXktO1Pa3youlzNziezS6jRXBgCQCAPIIGOMvNbdct96VYofexbBY7IsOZPnDbYG2CyjAQCZQhhARpiBPsXfelmmvSml/a3S6sHWgOLyNFcGAHgnwgDSyhgjb/92uTtek9x48gewncG+AXWz6BsAAGOEMIC0MX3dim9bJ9N5IKX9rYrxCs04TVZhSZorAwAcD2EAo2Y8T17zlsGphD03+QM4YTnTFskeN5XWAADIAsIAUmaMkWnbq/jODSlNHiRJVvUkhaafKitSlObqAACJIgwgJV5ni9wd62W6jz7xzwmFCxSavlh2DTMFAkC2EQaQFNPXqfiODSmPEpAke9xUOVMXyQpH0lgZACBVhAEkxET75O5ulLdvu5JZVGiEgmKFZpwmu3JCOksDAIwSYQDHZdyY3L2bB1cWTKVz4CH2xFlyGhbIcviRAwC/4Tczjsp4nrz9b8nd3SjFUpg98BCrqGxw8qCymjRWBwBIJ8IARkjHCAFJkmXJrp8rp/4kWbaTvgIBAGlHGMCwUY8QOMSqmqjQlIWymEoYAAKBMIBDIwRel2nfO6rjWKVVcqaeLLt8XJoqAwCMBcJAHjPRfrm73xjdCAFJKiwZbAmormcGQQAIIMJAHkrXCAGFCuRMnit7wgyWGAaAACMM5BETG5DXvE1u81YpnvoIAdmO7Imz5UyaIysUTl+BAICsIAzkAdPXJbdpi7wDO0bXEiDJHj9NTsN81hIAgBxCGMhhXler3L2bZNpG1zFQYoQAAOQywkCOGZonwN27adRDBCVGCABAPiAM5AjjxuUd2CF372ZpoGf0B2SEAADkDcJAwJlov9zmrfL2bZPi0dEfkBECAJB3CAMBZXo75TZtlndgp2S80R+QEQIAkLdyKgwYY2Q8T8YMTqBjJFmSZFmyDvsTVMYYmc4WuU2bZNqb03NQy5I9biojBADgHYwxw39kzNvXFGnwemLbgb6mHC7QYcAYIzcel+d58rzEPh1bti3bsmQ7juyAvJHGePJa98jbu1mmpz09B3VCsifMkFM3U1ZBcXqOCQABZowZvJ64rrzDPlieiG3bsm1bTigUiGvK0Vgm0f+tTwy9WUMhYDQsy5LjOL59A40bk7d/u9ymLdJAb3oOGimWM3GW7PHTuB0AAHr7g6XrugkHgGMZCgVB+bA5JFBhwHVdxaJp6CR3FLbjKBwOZ+zN8zpbZJfXnnA7Y4xMT7u8/dvlteyS3HhaXt8qqZQ9aY7s6no6BgKABn/fxmIxee7oJmM7lnAkIscJxhLugQgDxhjFotFRtwQkIhNvntuyS+7mFxQ+7UJZBSVH3cbEovJadsrbv12m92DaXtuqrJMzabas8nGBSqkAkEmZ/HB5ONu2FY5EfP/71/dhwI3HFYvFxvQ10/nmeW17Fd/0nGSMnIYFcibPHX5usEPgAXn73pLXtjc9owIkybJl106RPWm2bGYMBIBhY/nh8nDhcFhOyL/d9HxbmTFG8Xhcbjw9zeTJ8DxP0YEBRQoKRhUIvIP7Fd+0VjqUt9wDO2TXnyRF+wYnCNq/PX19ASTJCcuumznYKTBSmL7jAkAOMMYoOjAw6n4BqYjFYvKMUcinfdR82TJgjFE8FpObofs4ySgoKEjpHrvX1ar4G08dsTCQVVYr09WSrvIGFRTLmTh7sFOg49t8BwBZYzxPAwOjWK01TRzHUSiD/dNS5cswEIvFstIicCwFhYVJvXFeT4firz8puZm9vWGVVg1OEsSUwQBwTMYYDfT3Z7uMYU4opHDYX6O5fPcx0nVdXwUBSYpGo4ok2IfA9HUp3vh0RoOAVTVxMASU1RACAOA4jDGKjkFHwWS48fjgEEQfjTTwVRgY6tjhN8bz5LquQifo/GH6uxV74ykploGmKMuWPX6qnImzZRWVpf/4AJCDXNeVGePOgomIRaOyk2x1ziTfhAG/BoEh8VhseJapozF93Yq98aQU7Uvr61rFFbInTJdd2yArFEnrsQEgl3mep/gYj0ZLRiwa9c2wQ9+EgaHpH/0sFo0edYSB6etS7PUnpVia7kk5Ydm1DbLHT5NdWpWeYwJAHvH7B0xJw1Mf+2HIYfYr0NvDCP1uaCrkw+/zpDMIWOW1gwGgup5RAQAwCsmsLZBN8XhctuNkvXXAF1ecoLxp0mDHj6EwYPo6FXv9qdEFgXDh4KqB46fSFwBA3jPGpOXC6LeO6MdytA+Z2eCLMBCUN03S8AqJVv+hPgKj6CzoTF0ku24mawUAwCGt//VvKp63WEULT0/5d2MyK9n6weEfMrMl61chM0ZvWldXl5avWKEzzzxTS08/XXfeeWfKx4p3d4w6CEiSifUTBADgMPH9TWr9r39X8z//g3pfW5vSSICxnLDu41dcoYmTJumTV16Z8jE8z8v6iIesTzo0VmsPuK6rgYEBFRcXq7e3V0uWLtXTTz2lmpqapI9lSYrYg0MJD/+jvkN/J7rGQLhQ4SXvz/q9IgDwi+Z//o5i+3YP/zs0bqIq3n1pUi0FA/39Y3brefXq1eru7tavf/Mb3XfvvSkfJ9trF2T9NsFYNeU4jqPi4mJJUn9//6jWrTaSFCmUXVAkVYwb+ZwxUrTvUEDokYn2SQO9MtFemYFeaaDv7bAQ65c5uF9W5YRR/M8AIHfFDwy2FIT++j8JhQJjzJj2QTv33HP15JNPjvo4nucpmzcKAh8GPM/T4tNO0yWXXKLvffe7w48/9thj+shHP6q77rpLH7n8cklSR0eH3nfhhdqyZYu+/73vqba2NuXXbf/T/erf8GJK+1rhsJyCAtmFEXmvrFe8syvlOgAgl7g9nUd9PNFQkI4gkMx1JV2y3cchq7cJ0jVf9K9//Wtdd/312vjmm6qqqtJrr72m915wgb75zW/q61/72hHb79u3T5/45Cd1/333acKE1D6Vd//xN+p75dlRVg4ASMWxbh/E4/G0TDSUzHXlySef1M//9V9HdZtASn4dnHTKid5rn/jEJ1RTU6M77rhDu3fv1ocvv1xXfvKTRw0CkjRhwgSdvHChnl6zZmwLBQCkxVBLwTs7Gqbr822y15Wgy2oYSNebFgqFdP111+mOn/9cH778ci0+9VTdeuutI7bZt2+fOjsHm586Ozv19Jo1mjN7dlpeHwCQHUeEgjSNJEjkupJu2ezPn9XbBJ7nKZqm9aW7u7vVMGWKZs6cqdWrVqmkpGTE8+vWrdOXvvQlGQ1+wz//uc/pC1/4Quqv94ffqO9VbhMAgF9EGmao7OJPyB43KS3HO9F1RZIu/eAH9corr6inp0dVVVX67f33a+nSpSm9XqSg4Jjr32RaVjsQpvPeyNevu06S1NrSctTJG0477TStXbs2ba8nRgMCgC9EGmao/PwPqXDWAsXj8bRNZHei64ok/c/vf5+W15LSe01MVk70Gbjpppv06KOPavWqVYq7ru6+++5slwQAyLBIwwzVrvy6xn/h/6po9kJZlpW2z2n5dl3JahhIRwq666679NPbb9eDDz6oRYsW6SvXXqt/+slPxmAiI5oGACAVRlJTwylqHT9T3WXjFQsXKpn71UcLAUPSMatrtq4r2WwZyPoMhNGBgZTHV/75z3/Wx6+4Qvfcfbcuu+wySYOdA0+aO1c//MEPtHLlyjRWOpLV0iT3YFvGjg8AuWxNrEJ7vcLhfzvyVGR5irQ2KdLdpoK+zsE//YN/F/V0qKRu4vDtgGNdOEc7ZD1b1xXbthUpKMjIsROR9TAQj8VSWr543bp1et+FF+rGG2/UV669dsRzN998sx548EG98vLLGVv8obCoKCPHBYB8sLOjT3/d2pLUPhHHUllBWOUFIZUVhFRWGBr+ujBkDweE/r6+lGrK5nUlFAopFA6n/biJynoY8FxX0Wg0myUkLdsJDgCCzjNGD7y2V/3x9My8F7Kt4WCwoLZIZZFgdYmLRCKys7hyYda/W0FctS+bi0kAQC6wLUszao4cqpequGfU1hdTZVFY1aWFJ97BZ7J9Lcz6ldiyLIUCdnHN1jhQAMgVUddL+wXolInlOnVSReB+R4dCoayvXuuLq7ATCqXUbyAb/PCmAUAQGWPU0hvV5gM9equ9V3EvfXeph4KA9PaHzKBcV/zQ2pz9CjT4xjmOIzdN00hmkh/eNAAIkqjr6a22Xm060K22vvQPzzs8CAwJyodMx3F88QHTN1e2UDjs+zAQCod98aYBQBC09ES1qaVbb7WltxXgcEcLAtKh1oFwOC0rGGZSNkcQHM43YcCyLIUjEcV8OrLAsu2MDVMEgFwRO9QKsLGlW229mb0QnzqxXKccJQgMGWpxNinOZZNp4UjENx8wfRMGpENvnG2nPAlRJkVoFQCAY+qNuWrc16WNB7oVy1ArwOFOFASkwQ+ZkXBYA2laEC+dbJ99wMz6PAPvZIxRdGAgq0s5vlM4HKavAAAcxcH+mF7f16WtrT0agwwgKbEgcDg3Hh+DKeoTZ1mWIgUFvvqA6bsrnGVZikQivklyIYIAAByhpWdAG5q7tKMjtdn+UrVgQpkWTSxPah8nFJKRfNN/IOKj2wNDfNcyMMTzPMWi0ay2EITC4cDNgQAAmWKM0d7Ofm3Y16XmrrH/wDarpkTLplalfCGNx+NZDQRDfeP8OA+Cb8OAlN1bBuFIxFf3cwAgWzxjtL29V683d2VkaGAiplQW6dwZNbJH+Ynadd2sdFT3462Bw/k6DEiDgcCNx8dsvKhl2wqHw75MbgAwluKepy0tPXp9X5e6o5kb+m1Jqq8oVHlBSG/s7z7i+bqyAr131jg5dnoupJ7nKRaLjdkog1AoJMfnE9b5vg18aKyo7TgZf/NC4bBvJoAAgGwZiLt680C3Gvd3ayBNCwkdTXHY0ezaEs2uLVFJJKS452lza49i7tufUWuKIzp/Zm3agoB0aLG5SESu62b0tkGQPlz6PgwMGXrzPM9TPB5PaygIQmoDgEzricb1xr4ubWrpydgkQZJUX16oOeNKNbmicESzf8i2Nb2qWJtaeiRJFYUhvXd2rcJO+i+mQ1MWO46T9tZny7YVCoVk23ZgriuBCQPS29MWO44jz/Pkuq7cFN9A27blBOzNAoBM6OiLacO+Tm1r7VWmIkDRUCtATYlKC4596ZlVW6JNLT0qiTi6YPY4FYYy23drqPXZCYUGryvxeMpz3TiHwkUQWgLeKVBh4HC2bcs+1ARjjJHneTKeN/j3O7tBWJZsyxq88B+6+BMAAOS7/d0DWt/cqd0H+zP2GhPLCjR3fKkmVxQl1PmvtjiiurICnTWlSiWRsbtEHf5h0xgz+OfQNcUzRnrHdcU67JqSCx8qfd+BEACQXi09A3px90Ht687M8EBL0rSqYi2sK1N1cSTp/V3PpLWPAE6MMAAAeaJ7IK51ew/qrbbejBzfsSzNqi3RggllKjvOrQD4D+8WAOS4qOtpfVOn3tjflZEpgyOOrbnjSzVvXKkKw8zPEkSEAQDIUZ4x2nSgW680dWZkiGBJ2NH8CWWaXVuSkR7/GDuEAQDIMcYY7TrYr5d2d6hzIP0TtlUUhnRyXbmmVxePekZA+ANhAABySGtPVC/u7lBzBjoHji+JaGFduSZXFAa+9zxGIgwAQA7oica1bs9BbctA58DJFYU6ua5c40sL0n5s+ANhAAACLOp62tDcqTf2dctN4+AwS9KMmmItmFCuqqJw2o4LfyIMAEAAecZoc0uPXtl7UP1p7BwYsi3NqS3R/AllYzrpD7KLdxoAAsQYoz2d/Xpxd4cO9qevc2DYtrRgQpnmji9VQYanAIb/EAYAICDaegc7BzZ1pa9zoCVpzrgSnTKxQkXMEZC3CAMA4HO9UVfr9nZoa2t6OwdOrijUkvpKVdInIO8RBgDAp4wx2tLaoxd2dSiWxqkDq4vCWjq5UhPLC9N2TAQbYQAAfKh7IK5ndrSl9ZZAcdjR4voKzawuZp4AjEAYAAAfMcZoY0uPXtrdoXiaWgNCtqWFdWVaMKFMIZtpg3EkwgAA+ETnQFzPbm9L2+yBlqRZtSVaPInOgTg+wgAAZJlnjN7c362X9x5MW2tAfXmhlkyuUFVRJC3HQ24jDABAFh3sj2nN9jYd6Imm5XhVRWEtmVypejoHIgmEAQDIAs8YvbGvS6/s7UzLNMJFIXuwc2BNCSsJImmEAQAYY+19g60Brb2jbw0IHZo5cMGEMoUdOgciNYQBABgjnjHa0NylV5sOKh1dA2bVlOjUSeWsIYBR4ycIAMZAW29Ua7a3qa0vNupjlRWEtGxqlerK6BeA9CAMAEAGuZ7Ra82dWt/UqXSME5g/vlSL6yuYLwBpRRgAgBT09PSoo6NDZWVlKi8vP+o2LT2DrQEd/aNvDSgvCGn5tGqNLy0Y9bGAdyJaAkCCtm3bpmuuuUbz589XWVmZJk+erHHjxumee+4ZsZ3rGb20u0N/enPfqIOAJWnhhDJ9cH4dQQAZYxmThjEtAJAHLr74Yr3++uu69NJLtXTpUk2cOFH/+I//KGOMHn/8cUlS10Bcq7a1qK139K0BlYWDrQG1JYQAZBYtAwCQoL179+rSSy/Vl7/8ZTmOo4suukiTJk1SPB6XJO0+2Kc/NDaPOghYkhZNLNcl8+oIAhgThAEASND8+fN1xx13aOHChVq5cqU8zxt+bt2eDj2+pUVRd3SNrdVFYV0yb4IWT6qQYzN5EMYGYQAAEvSLX/xC9913n772ta+NeLy9L6b1zV2jOrZtSadOKtcH5k1QdTHrCWBsMZoAABJUWlqqT3ziE+rr6xvxeNT1jrFHYmqKI1o+rYpFhZA1hAEASFK6+l3blrR4UoXmTyhjPQFkFWEAAJIQcz1tauke9XHGlUS0fFq1KgrDaagKGB3CAAAkqKMvpj+99pa27Woa8Xi0v1/tB5pVNa7uhMdwLEun1Vdo7vhSWgPgG4QBAEjAW229+rf7f6efXPc5RQf6Rz73xqv6+sVn6F2XfULXfPuWYx6jtjiic6ZXq5zWAPgMkw4BwHG4ntGLuzv05oFu3fndb6h183rdfvvtkqTzzz9fjY2Nampq0h//+Ef99Ke361drtx31OHPHlWrp5EqGC8KXaBkAgGPoica1elurDvREJUnFpeXa2NKiBx54QJKG/5akxsZGFZWWHXGMkG3p7KlVmlFdMjZFAykgDADAUezt7NeTb7VqIP72sMH3fHylDuzZpUdXrzli+3CkQJ/55vdGPFZeENJ5M2tVVcRtAfgbtwkA4DDGDC45/MrezlEdZ2pVkZZPrVbYYW43+B8tAwBwyEDc1VNvtWlPZ/+JNz4GS9LSyZWaN75UFqMFEBCEAQDQ4JTCT2w5oO6om/IxisKOzptRw1LDCBzCAIC8t69rQE9sPTCqRYbqygr0ruk1Kgo7aawMGBuEAQB5bUd7r558q1XeKHpPnVxXplMnVTCJEAKLMAAgb715oFtrd7anvH/YsXTOtBo1VBalsSpg7BEGAOQdY4xe2dup15pTHzFQXRTWuTNrVV7Ar1EEHz/FAPKKZ4ye3dGuLa09KR9jdk2JzphSqZDNsEHkBsIAgLwR9zyt3taq3QdTGzpoW9JZU6o0u7Y0zZUB2UUYAJAX+uOuntjSMjy1cLJKI47Om1mrmuJImisDso8wACDndQ/E9djmA+ociKe0/6TyQr1reo0KQtwWQG4iDADIaW29Uf3vlgPqi3kn3vgoZtYUa9nUaoYNIqcRBgDkrKaufv11S4tiKU4icHJdmRZPqmBaYeQ8wgCAnLS9rVdPbU99MqEzGio1b/yRSxIDuYgwACDnNO7v0vO7OlLa17akc6bXaFpVcXqLAnyMMAAgZxhjtG7PQW3Y15XS/mHb0vmzalVXVpjmygB/IwwAyAmeMXpme5u2tvWmtH9R2NF7Z9WqmqGDyEOEAQCBF3M9rdrWqr2dqU0mVFEY0ntnjVMpUwsjT/GTDyDQoq6nxzYdUEtvapMJjSuJ6PxZtSoMsfQw8hdhAEBgxVxP/7s59SAwuaJQ586oYY0B5D3CAIBAinveqKYXnl1borOmVDGZECDCAIAAcj2jVVtb1dw9kNL+p0ws1ykTy5lMCDiEMAAgUDxjtPqtVu1JobOgJenMKVU6aRyrDgKHIwwACAzPGD31Vpt2dfQlva9tSe+aXqOpTCYEHIEwACAQjDF6ZkebtrcnP49AxLF0/qxxmlBakIHKgOAjDADwPWOM1u5q19bW5INAcdjRe2ePU1VROAOVAbmBMADA14wxenF3hzYe6El63+Kwo4tOGq8yJhMCjovBtQB87ZWmTr2xvzvp/QpDtt43ZxxBAEgAYQCAb61v6tRrTZ1J71fgDAaBikJuDQCJIAwA8KU39nVp3d6DSe8XdixdMGecqopYcAhIFGEAgO9sPNCtF3Z3JL1fyLb03lnjVMPKg0BSCAMAfGVra4+e29me9H6OZek9s2o1nuGDQNIIAwB8Y3t7r9Zsb0t6P9uS3j2zRnVlhRmoCsh9hAEAvrCzo09PbmuVSXI/S9K5M2pVX1GUibKAvEAYAJB1ezr7tXpbS0pB4JzpNZpSSRAARoMwACCrmrv69dctLfKSTQKSlk2r1vRq1hoARoswACBrOvtj+uvWVrkm+SRw1pQqzaopyUBVQP4hDADIimjc0+NbWhR1vaT3XTq5kmWIgTQiDAAYc54xWv1WqzoH4knvu3hShRZMKMtAVUD+IgwAGHMv7OrQ3s7+pPc7ua5ciyaWZ6AiIL8RBgCMqY0HuvXmgeQXHpo/vlSLJxEEgEwgDAAYM02d/VqbwuyCc2pLtHRypSzLykBVAAgDAMZEZ39Mq1KYVGhmdbHOmlJFEAAyiDAAIONSHTlQV1agZdOqCQJAhhEGAGRUqiMHygpCOm9GjWyCAJBxhAEAGZXKyIGwM7gCYUHIyVBVAA5HGACQMamMHBhaeKiiMJyZogAcgTAAICNSHTlwekOl6stZihgYS4QBAGmX6siBObUlmss0w8CYIwwASKvRjBw4kyGEQFYQBgCkzWhGDpzLyAEgawgDANJmNCMHChk5AGQNYQBAWjByAAguwgCAUWPkABBshAEAo9I1EGfkABBwhAEAKXM9o1XbUhg5UMrIAcBPCAMAUvbi7g619caS2qesIKRzZzJyAPATwgCAlOxo7026w2DYtnT+TEYOAH5DGACQtK6BuNZsb0tqn8GRAzWqLGLkAOA3hAEASXE9o9XbWhTzkusyuHRypeorijJUFYDRIAwASMpLezrUmmQ/gTm1JZo3npEDgF8RBgAkbEd7rxr3J9dPoLY4ojMaGDkA+BlhAEBCugbiWrMjuX4CEcfSuTNq5NgEAcDPQtkuAMDYM25c3sEWKdon47mSJMt2pEiR7IpaWc7IXw2D/QRaFXOT6yewfFqNSgv4NQP4HWcpkAe83i65u9+U19Yst3WPTFerZI5xYbcsWWU1cmrqZVfXyZk8V7v6HbX2RpN6zfnjSzWlkg6DQBBYxhzrNwKAIDOeJ7d5m+JbXpK7Z8vgg5YlmQRnC7Ts4cDg1M/W3vqlWtsVUTSB1oHa4oguOmk8tweAgCAMADnGGKP49vWKvfqETF/3oYt6ctMFH+HQMfpKx+n5KRfoQPzYcwVEHEuXzKtTGbcHgMAgDAA5xOvu0MDzf5C3b3vmXkOWNk55l94onHbUxYnePbNGUyqLM/b6ANKPMADkAGOM4pteUPTVJyTPPXZ/gHSxLLUU1+n5ye9Rr3l7auF540t1RkNVZl8bQNoRBoCAM56ngRf+KHfbq2P+2lE7opdnXKxddqVqiiO6mH4CQCARBoAAM56rgacflrtnY/ZqkLRz6gpNXny2yosKslYHgNQRBoCAMp6ngWd+J3dXY7ZLkSQ5U+ar4OzLZNnMZQYEDWctEFCxxmd8EwQkyd35hmKNz2S7DAApIAwAAeS2NSm2fnW2yzhCbP1quW3N2S4DQJIIA0DAGDeugWd+l+0yjmngmYdl3Hi2ywCQBMIAEDCx9atlutozP3wwFcbIdLUrtv7JbFcCIAmEASBATLRfsY3PS0ed7scvjGKbnpeJ9me7EAAJIgwAARLf9urgpEJ+58YHawUQCIQBICCM5ym2cW22y0hYbONaGW+UayIAGBOEASAg3OZtMr2dY/Z6vQNRzfv//lH/9zePprS/6e2U27wtzVUByATCABAQbvNbg6sHjpFb/nu1ls6cnPoBLHuwZgC+RxgAAsJr2T36pYgTtKWpVZv2tujCU+ekfhDjyWvdk76iAGQMYQAIAON58tpHN5mP53lafP1P9ff3/XnE4//76mZVffpGPfzchuHHvnXvo7rpExeM6vUkyWtvpt8AEACEASAATGfLqEcR2LatGz70Lv3yf19Qe3efJGn9jiZ9+vbf6sYr3qvLz1ooSfrDi42aVVej2RNrR1233Phg7QB8LZTtAgCcmNfdkZbjXLF8kX7w0F/18z8/q5XnLdFHfvxrfWLFKfrbS1YMb/PCll168Nn1+t3a19XdH1XcdVVWVKBvXv7u1Grv6ZBdOT4t9QPIDFYtBAIgvvMNDax5OC3H+tXjL+im3z6miVXlmjquUvddd6WcY6w0+OvV6/TG7v36/qcuSvn1CpZfrtCU+SnvDyDzuE0ABEEa77tfsXyRegdiMjK66ysfP2YQSJsgTJIE5DluEwBB4DhpO9T1d/9RktTa1SvHto677VXnnjb6F3TCoz8GgIyiZQAIglAkLYe5+b/+V4++vFF/vfkLirue7lm1Li3HPR4rTbUDyBzCABAAduWEUR/j7r++qH/+0zN64IardPLUibr24rN12x+eUiye2WZ8Og8C/kcYAALALiqVVViS8v5/eWWTrrvrD/rllz+iM2Y3SJK++L6z1NU7oPuefiVNVR7JKiyRVVSaseMDSA/CABAQdk29pOPf4z+al7ft0ad/+lv9v09eqA+dsWD48fLiQv2fC8/SP/3+KbmZmBjIsmTXjmI6YwBjhqGFQEBEX1+j2PpVUlBOWctSeNF5isxfnu1KAJwALQNAQISmzAtOEJAkYxRqmJftKgAkgDAABIRdVi174kzJSv5WwZizbNkTZ8ouq852JQASQBgAAiR80hnBaB0w3mCtAAKBMAAEiFM3Q1ZplVLpSDh2LFmlVXLqZmS7EAAJIgwAAWJZliJLLpTk59YBo8iSi2QF4XYGAEmEASBwQpNmKTRzsfzZOmApNOs0hSbNzHYhAJJAGAACKLL4AlnF5f7qTGhZskrKFVn83mxXAiBJhAEggKxwRAXLPyx/tQ5YKlj2YdYiAAKIMAAElFM7WQUrPuqP1gHLUsE5H5XDjINAIBEGgAALTZ6jguUfkSxb2WklsCTLVsGKjypUPycLrw8gHZiOGMgB7v4d6l/9W8mNjd08BJYlOWEVnnuFnPFTx+Y1AWQEYQDIEV5vl6IvPiJ3z6YxeT2nfo4iSy+WXVw2Jq8HIHMIA0AOMcbI3fWmBl74kxTtV/rnI7CkSKEKzng/6w4AOYQwAOQgM9Cn2OYXFd/8kkx/92CTfqqn+qF9rcJShWYvUXjOUlmRovQWDCCrCANADjOeJ7dpq+Jb1sndu0XDLQWWLRnv6DuNeM6SM2mWQrOXDE6FbNPnGMhFhAEgT5h4VF7HfnntzfLamuW1NclE+2TcuCTJckKyIkWyqyfKrq6TXVUnu3I88wYAeYAwAABAnqPNDwCAPEcYAAAgzxEGAADIc4QBAADyHGEAAIA8RxgAACDPEQYAAMhzhAEAAPIcYQAAgDxHGAAAIM8RBgAAyHOEAQAA8hxhAACAPEcYAAAgzxEGAADIc6FsF3C4yOJrZIcismxHlu3ICb/9tWXbbz/nOLJDEdnDzzlHPGfZjmzbkmVbchxb1ju+tm1LtmMNb3Pc5yxLTsiWY1tybEuRQ1+Hhv/tvP2c8/Z2ocO2dY72tWXJtiw5lhR27OGvQ44tx9Lgv21LYds6yteDz4dte/hrx7JkWZJtSZalQ8eXLEmObcmWBv8vtoa/ti3JsQ7/evAYljGS8WR5cWnE197gH+/Yz1nGk1z37a+9uOS5Mp4nxaMyrit53uBj8ZiM5w5+HYtJQ18PbTu0XSz69j6eKy8Wl3E9Gc+TF43Lcwf3Ma4nLxaX5779tTn0tRuLyxy2nRuNH/a1K+MZea459O9D+3tm8DnXyLhGnuvJjXmHjmnkxtxD+7y9n2eMXGMU9Yxco3d8/c5/D37tafBr1+jQc29//a9me1bPy3Th/Ob85vz27/lNywAAAHmOMAAAQJ4jDAAAkOcIAwAA5DnCAAAAeY4wAABAniMMAACQ5wgDAADkOcIAAAB5jjAAAECeIwwAAJDnCAMAAOQ5wgAAAHmOMAAAQJ4jDAAAkOcIAwAA5DnCAAAAeY4wAABAniMMAACQ5wgDAADkOcIAAAB5jjAAAECeIwwAAJDnCAMAAOQ5wgAAAPnO5Kj+/n7zne98x/T392e7lCP4uTZjqG80/FxbLvHz99nPtRlDfaPh59pGyzLGmGwHkkzo7OxURUWFDh48qPLy8myXM4Kfa5OobzT8XFsu8fP32c+1SdQ3Gn6ubbS4TQAAQJ4jDAAAkOcIAwAA5LmcDQMFBQX6zne+o4KCgmyXcgQ/1yZR32j4ubZc4ufvs59rk6hvNPxc22jlbAdCAACQmJxtGQAAAIkhDAAAkOcIAwAA5LmcCwM33HCDzjnnHH3qU59SNBod8VxfX58uueQSnXvuubrgggvU1tbmq/qG/OAHP9DSpUuzXlM8HtfVV1+tc845R3/7t387ZvUkWt+Qsf5+He5YtfnhZy0XcX6nrybO7xPLp/M7p8LAyy+/rObmZj311FOaP3++HnzwwRHPP/LII1q4cKFWr16tj3/84/rP//xPX9UnSV1dXdqwYYMvavqf//kfTZ48WU899ZR6e3v1zDPPjFldidQnjf33K9Hasv2zlos4v9NbE+d36rVl+2ctE3IqDDz77LN63/veJ0m66KKLjvjhnj17tnp7eyVJHR0dGjdunK/qk6Sf/vSnuvbaa31RUyL1ZrM+aey/X4c7Xm3Z/lnLRZzf6a2J8/v48u38DmW7gHTq6OjQpEmTJEkVFRVHNN3MnDlTGzZs0MKFC2VZltauXeur+g4ePKj169fr29/+ti9q6ujoGJ5/+2j1Zru+bHy/Eq0t2z9ruYjzO701cX6nXlu2f9YyIZAtA83NzVqxYsURf4wx6uzslDT4RlZXV4/Y75577tF5552nDRs26KabbtLNN9/sq/puu+02feUrX8lITcdSVVV1zJqO95wf6svG9+twx6ttrH7WchHnd/pwfqcu387vQIaBuro6Pf3000f8ef/736+//OUvkqQ///nPWr58+RH7Dr2hlZWV6ujo8FV9W7Zs0fe+9z1ddNFF2rx5s374wx9mpL7DnXXWWces6XjPjZXj1ZCN71eitUlj87OWizi/04fzOzO1STl4fmdv9eTMuP76682KFSvMlVdeaQYGBowxxnzhC18wxhhz8OBB8/73v9+ce+65Zvny5Wbjxo2+qu9wS5YsyVpNQ/XEYjHzmc98xqxYscJ89atfHbN6Eq3vcGP5/TrcsWrzw89aLuL8Hn1NnN+Jy6fzm+mIAQDIc4G8TQAAANKHMAAAQJ4jDAAAkOcIAwAA5DnCQB64++67VVlZmZZjbd++XZZlKRQKac+ePSOea2pqUigUkmVZ2r59+4jnHnroIZ133nmqqKhQaWmpFi1apJtvvnl4Io901gjkm6uvvlqWZemLX/ziEc99+ctflmVZuvrqq4cfa25u1le/+lXNmDFDBQUFamho0KWXXqrHH398eJtp06bptttuG4Pq4QeEAaRk0qRJ+o//+I8Rj91zzz2qr68/YttvfetbuuKKK3T66afrkUce0YYNG3Trrbfq1VdfzYk5vQE/aGho0P3336++vr7hx/r7+3XfffdpypQpw49t375dS5Ys0RNPPKFbbrlF69ev16OPPqp3v/vdWZv6F9lHGAiARx99VCtWrFBlZaVqamp0ySWXaOvWrZKkVatWybKsEZNevPLKK8OfzletWqXPfvazOnjwoCzLkmVZuvHGGyVJ7e3t+sxnPqOqqioVFxfr4osv1ubNmxOqaeXKlbrrrrtGPHb33Xdr5cqVIx57/vnn9f3vf1+33nqrfvzjH2vZsmWaNm2aLrjgAj300ENHbA8gNaeddpqmTJmihx9+ePixhx9+WA0NDVq8ePHwY0MtBc8//7w++tGPas6cOVqwYIGuu+46Pffcc9koHT5AGAiAnp4eXXfddXrhhRf0+OOPy7ZtffjDH5bneSfcd9myZbrttttUXl6upqYmNTU16YYbbpA02LT44osv6ve//72effZZGWP0/ve/X7FY7ITH/eAHP6j29nY9/fTTkqSnn35abW1tuvTSS0ds95vf/EalpaX68pe/fNTjcGsASJ/PfvazI0L6nXfeqWuuuWb4321tbXr00Ud17bXXqqSk5Ij9OR/zV04tVJSrPvKRj4z4969+9SuNHz9eb7zxxgn3jUQiqqiokGVZqqurG3588+bN+v3vf681a9Zo2bJlkgYv3A0NDfrv//5vfexjHzvuccPhsK666irdeeedWrFihe68805dddVVCofDI7bbvHmzZsyYccTjANLv05/+tL75zW8O9+1Zs2aN7r//fq1atUrS4BS/xhjNnTs3u4XCd2gZCICtW7fqyiuv1IwZM1ReXq7p06dLknbu3JnyMRsbGxUKhXTmmWcOP1ZTU6OTTjpJjY2NkqSLL75YpaWlKi0t1YIFC444xt/8zd/ogQceUHNzsx544IERn0CGGGNkWVbKdQJIXG1trT7wgQ/onnvu0V133aUPfOADqq2tHX5+aMJZzkm8Ey0DAXDppZeqoaFBv/jFLzRp0iR5nqeFCxcqGo2qtLRU0tsnuaSEmvmPNQv14RfvX/7yl8OdkY72yX7hwoWaO3euPvnJT2revHlauHChXnnllRHbzJkzR08//bRisRitA8AYuOaaa4ZX+/vZz3424rnZs2fLsiw1Njbqsssuy0J18CtaBnyutbVVjY2N+va3v633vOc9mjdvntrb24efHzdunKTBYX1D3nlBjkQicl13xGPz589XPB4fsQ53a2urNm3apHnz5kmS6uvrNWvWLM2aNUtTp049an3XXHONVq1addRWAUm68sor1d3drTvuuOOoz+fEal+Aj1x00UWKRqOKRqO68MILRzxXXV2tCy+8UD/72c/U09NzxL6cj/mLMOBzVVVVqqmp0b//+79ry5YteuKJJ3TdddcNPz9r1iw1NDToxhtv1KZNm/THP/5Rt95664hjTJs2Td3d3Xr88cfV0tKi3t5ezZ49Wx/60If0+c9/Xk8//bReffVVXXXVVaqvr9eHPvShhOv7/Oc/rwMHDuhzn/vcUZ8/88wz9Y1vfEPXX3+9vvGNb+jZZ5/Vjh079Pjjj+tjH/uY7rnnntS+MQCOynEcNTY2qrGxUY7jHPH8HXfcIdd1dcYZZ+ihhx7S5s2b1djYqNtvv11nn312FiqGHxAGfM62bd1///166aWXtHDhQn3961/Xj3/84+Hnw+Gw7rvvPr355ps65ZRT9KMf/Ujf/e53Rxxj2bJl+uIXv6grrrhC48aN0y233CJJuuuuu7RkyRJdcsklOvvss2WM0Z/+9KekmvNDoZBqa2sVCh37jtOPfvQj3XvvvVq7dq0uvPDC4WFMixYtYmghkAHl5eUqLy8/6nPTp0/XunXr9O53v1vXX3+9Fi5cqAsuuECPP/64fv7zn49xpfALljAGACDP0TIAAECeIwwAAJDnCAMAAOQ5wgAAAHmOMAAAQJ4jDAAAkOcIAwAA5DnCAAAAeY4wAABAniMMAACQ5wgDAADkuf8fOQ13pgiuoUUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pc_alpha_bootstrap = [0.001, 0.01, 0.05, 0.1, 0.2] # This can be adapted \n",
    "boot_samples = 200\n",
    "\n",
    "# The block-length of the bootstrap can optionally be used to better deal with autocorrelation, \n",
    "# but its effect was not yet evaluated.\n",
    "boot_blocklength = 1\n",
    "\n",
    "## Create PCMCI object to call run_bootstrap_of\n",
    "pcmci = PCMCI(dataframe=dataframe,\n",
    "        cond_ind_test=ParCorr(),\n",
    "        verbosity=0,\n",
    "        )\n",
    "\n",
    "# Call bootstrap for the chosen method (here 'run_pcmciplus') and pass method arguments  \n",
    "results = pcmci.run_bootstrap_of(\n",
    "        method='run_pcmciplus', \n",
    "        method_args={'tau_max':tau_max, 'pc_alpha':pc_alpha_bootstrap}, \n",
    "        boot_samples=boot_samples,\n",
    "        boot_blocklength=boot_blocklength,\n",
    "        seed=123)\n",
    "\n",
    "# Output graph, link frequencies (confidence measure), and mean test statistic values (val_mat)\n",
    "boot_linkfreq = results['summary_results']['link_frequency']\n",
    "boot_graph = results['summary_results']['most_frequent_links']\n",
    "val_mat = results['summary_results']['val_matrix_mean']\n",
    "\n",
    "# Plot\n",
    "tp.plot_graph(\n",
    "    graph = boot_graph,\n",
    "    val_matrix= val_mat,\n",
    "    link_width = boot_linkfreq,\n",
    "    var_names=dataframe.var_names,\n",
    "    ); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrap aggregation with optimized ``pc_alpha`` for CMIknn\n",
    "\n",
    "Finally, we showcase the Bagged-PCMCI+ approach with optimized ``pc_alpha`` and the nonlinear conditional independence test ``CMIknn``. See the corresponding [tutorial](https://github.com/jakobrunge/tigramite/blob/master/tutorials/causal_discovery/tigramite_tutorial_conditional_independence_tests.ipynb). Using the standard initialization ``CMIknn(significance='shuffle_test')`` would result in each test carrying out a computationally intensive permutation-testing scheme to obtain the null distribution. An alternative with less statistical rigour is to just make test decisions based on a fixed threshold on the conditional mutual information using ``CMIknn(significance='fixed_thres')``, see the explanation in the tutorial. This is illustrated here.\n",
    "\n",
    "We create nonlinear data as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the time series length\n",
    "T = 500\n",
    "\n",
    "# Specify the model (note that here, unlike in the typed equations, variables\n",
    "# are indexed starting from 0)\n",
    "def lin(x): return x\n",
    "def nonlin(x): return .2 * (x + 5. * x**2 * np.exp(-x**2 / 20.))\n",
    "\n",
    "links = {0: [((0, -1), 0.3, lin), ((2, 0), 0.5, lin), ((3, -1), -0.7, nonlin)],   # X1\n",
    "         1: [((1, -1), 0.3, lin)],                                             # X2\n",
    "         2: [((2, -1), 0.3, lin), ((1, -2), 0.4, nonlin)],                        # X3\n",
    "         3: [((3, -1), 0.3, lin)]                                              # X4                            \n",
    "        }\n",
    "\n",
    "# Generate data according to the full structural causal process\n",
    "data, nonstationarity_indicator = toys.structural_causal_process(\n",
    "    links=links, T=500, seed=seed)\n",
    "# Initialize dataframe object, specify variable names\n",
    "dataframe = pp.DataFrame(data, var_names=var_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a range of fixed thresholds, these are used as pc_alpha (with a slight abuse of the parameter name)\n",
    "# This can be adapted, higher thresholds lead to stricter link decisions and, hence, sparser graphs\n",
    "fixed_thresholds = [0.01, 0.025, 0.05, 0.1]\n",
    "boot_samples = 100\n",
    "\n",
    "# The block-length of the bootstrap can optionally be used to better deal with autocorrelation, \n",
    "# but its effect was not yet evaluated.\n",
    "boot_blocklength = 1\n",
    "\n",
    "## Create PCMCI object to call run_bootstrap_of\n",
    "pcmci = PCMCI(dataframe=dataframe,\n",
    "        cond_ind_test=CMIknn(significance='fixed_thres', model_selection_folds=5),\n",
    "        verbosity=0,\n",
    "        )\n",
    "\n",
    "# Call bootstrap for the chosen method (here 'run_pcmciplus') and pass method arguments  \n",
    "results = pcmci.run_bootstrap_of(\n",
    "        method='run_pcmciplus', \n",
    "        method_args={'tau_max':tau_max, 'pc_alpha':fixed_thresholds}, \n",
    "        boot_samples=boot_samples,\n",
    "        boot_blocklength=boot_blocklength,\n",
    "        seed=123)\n",
    "\n",
    "# Output graph, link frequencies (confidence measure), and mean test statistic values (val_mat)\n",
    "boot_linkfreq = results['summary_results']['link_frequency']\n",
    "boot_graph = results['summary_results']['most_frequent_links']\n",
    "val_mat = results['summary_results']['val_matrix_mean']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAHBCAYAAAD0E7h1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6fUlEQVR4nO3deZwcdZ3/8Xf1PdNzZxJyQojhjigKggFccEVECKxAQA65EeQQQRZx8WB1QZEfq8ICitwosIug4nKpsBGCWUEQ1CVKCCYESEgyR+bsnu6u7++P6gyZzEym76qaej0fjzwy6a7j06nprndXfQ/LGGMEAAACK+R2AQAAwF2EAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgICLuF0AgNoyxpZyKSkzIJmcZGznCSskWWEpWi+FE7IsvisAQUEYACY5Y4yU6ZVSnVK2T8qmJJkJ1rJkIgkp0iAl2qRooyzLqkW5AFxgGWMm+lQA4EPGzkiDG6XB9ZI9VN7GQjGpbppU1y4rFK1MgQA8gzAATDLG5KT+t6WBddXZQf10KTlTlhWuzvYB1BxhAJhEzFCP1PP38q8ETCQUk5p2lBVrqu5+ANQEYQCYBIyxpd7VUmpjbXecaJcad6CxIeBzhAHA54ydlbpXOI0D3RBpkFp2khWiPTLgV4QBwMeMnZW6/irlBt0tJFwnte5KIAB8imt7gE8ZOyd1/839ICA5NXT/zakJgO8QBgC/6ntDyg64XcW7sgNOTQB8hzAA+JBJd9e+sWAhUhud2gD4CmEA8BljZ5zug17V83enRgC+QRgA/Kb3Dclk3a5ifCYr9a5xuwoARSAMAD5ickNSutPtMiaW7nBqBeALhAHATwbXu11B4fxUKxBwhAHAJ4yx/XWCHVzv1AzA8wgDgF+kuyRT/X78a956Rwd/6lztceDxet9BJ+qBh39T2oZMzqkZgOcxAiHgE6Z3lTMlsar7ll37zka9s6FT71+ws9Zv6NQHDzlFf332ASWTdcVvrG6qrMa5Fa8RQGUxdijgF5k+VTsISNKM7do1Y7t2SdK0qW1qa2lSZ3dPaWEg49J8CQCKwm0CwAeMsaVsecMO27at3fZfrC9984YRjz/xP8sUn71wzNsBf3jpFdm2rTmztittp9lB2g0APkAYAPygAsMOh0IhXf75U/WDOx9SV3ePJOnl/3tVx531L7rqX87T4iM/NmL5js5unXrhv+qH1/1LeTv20pDJAMZEGAD8oEIj+p10zCfU3tas62/9T7359js64qRLdPLiw3TpeSePWC6dHtLRp1+myz9/qhbus2d5O2U0QsDzaEAI+IBJdUg9r1dkWz+8+yFdcfVNmjV9mubOmaGH7vyOwuHwu/syRied+1XtPH97XfnPny1/h03vkZVoK387AKqGMAD4QCXDQF//gKbtfqjmz52tZY/ePqph4NLfv6R/OOoc7bn7/OHH7v6Pf9V7t/h3UZrmyUpMKadkAFVGbwLAF6yKbenCL18rSdrY2a1wePSdwgP2fb9y635fsf3J4m4k4HW8SwE/CFUmt3/12z/QI795VssevU3ZXE633ftwRba7TVZ44mUAuIowAPhBpL7sTdz645/r33/wEz1893V63x4766KzP61rb7xHmUyVZ0CsQO0AqoswAPiAFYpIoVjJ6z/25O90wZev1T03fkP77f1eSdKFZx2nnt5+3fPAo5Uqc7RQzKkdgKcRBgC/iDaUtNoLLy/XcWd9Wdd89UIdffjBw483NTbogjOP0zU33KVcrkpzHpRYM4DaojcB4BNm4B2p7w23yyhOw/ay6kscvRBAzXBlAPCLRJsq2aug+qx8zQC8jjAA+IQVikpxH51c421OzQA8jzAA+ImfLrn7qVYg4AgDgI9Y0aQUSbpdxsQiSadWAL5AGAD8pnGu2xVMzA81AhhGGAB8xorWS8lZbpcxvuRsp0YAvkEYAPyofoY3bxdEklL9dLerAFAkwgDgQ5ZlSU3zvDXuvxWWmuc5tQHwFcIA4FNWJCG17OqNWQGtkNSyq6xwwu1KAJTAA58iAEoWqZMS7ZLlYn9+K+IEAdoJAL7FDCKAn6U2SNl+KRKTTETKDtZ2/+GE1LITVwQAn2NuAsCnTKZv9FwFVlga6pVUg7d1/UwpOUOWF25TACgLYQDwIWNnpJ7XJTPWbIOWZGwpO1CdnUfqpaYdZUW4LQBMFoQBwGeMMVLf6olP9uE6KZeWMr2V2XG0SaqfJsVa6DEATDKEAcBnzOB6KbWxsIXDcal+lrN8qkOyh4rbWSgmJaZIdVNlhePFFwvAFwgDgI+M2U5gW+qny9pipkNj55wrCtl+p7GhnXVuKUhO98BQxOmhEElKkXpZIQ+NYwCgaggDgE9su53AGKJNUnIWl/QBTIhmwIAPGGOk/rcKDwKhWL6lP0EAwMQIA4AfpDYU0TvAyl8R4BI/gMIQBgCPM0M9hTcYlKT67WRF6qpXEIBJhzAAeJjJDjq3BwoVbZJirdUrCMCkRBgAPMrYWalvjQoeTTAUpZ0AgJIQBgAPMsZ2goDJFriGJSVn004AQEkIA4DHGGOkgbelXBGTDtXRTgBA6QgDgNekNkpDPYUvH22S4rQTAFA6wgDgIU7PgQ2FrxBOSMmZtBMAUBbCAOARRfccsCJSwxymEAZQNj5FAA8wdqa4ngOynCAQilazLAABQRgAXFZ8zwE5twZoMAigQggDgIucOQfelnKpwldKtMuKNVevKACBQxgA3JTaKGWK7DmQmFq9egAEEmEAcAk9BwB4BWEAcAE9BwB4CZ8sQI3RcwCA1xAGgBoqrefALHoOAKgqwgBQI6X1HJgqK9ZUvaIAQIQBoHZK6jnQXr16ACCPMADUgEl303MAgGcRBoAqM0ObnCmJC0XPAQA1xqcNUEVmqLe4LoT0HADgAsIAUCUm0yf1v1ncSvQcAOACwgBQBSbTX+RYAqLnAADXEAaACjPZAanvDRUVBOg5AMBFhAGggkx2UOotMgiE6+g5AMBVhAGgQkw2JfWtlmQXvlI4ITVsT88BAK7iEwioAJNLO0HAFBEEQnEnCITC1SsMAApAGADKZHJDUu9qyeQKXykUkxp3kBWKVK8wACgQYQAogxMEVhU38VAoShAA4CmEAaBEzlTEq4sLAlZEatiBQYUAeAphACiBsbPOrQE7U/hKVkRqnCsrHKteYQBQAsIAUKR3g8BQ4StZYefWAEEAgAcRBoAiGDvnDChkpwtfyQo5twbC8eoVBgBlIAwABTImHwRyqSLWygeBSKJqdQFAuQgDQAGMsZ25BnKDRaxlSY3bM/EQAM+jbxMwAeeKwBopO1DEWpYzoFCkvmp1AUClEAaAbTB2toRbA5bUMEdWNFm1ugCgkggDwDicIYbfKK77oCQlZ8uKNlSnKACoAsIAMAZnGuI1xQ0xLDlBINZYnaIAoEoIA8BWzFCv1P+mipqGWJLqZ8qKNVWlJgCoJsIAsAWT7pIG1ha/Yv0MWfGWitcDALVAGAAkGWOk1EYptaH4leumy4q3Vr4oAKgRwgACzxjjXA0Y6i5+5brpshJtFa8JAGqJMIBAM8Z22gdk+opc05KSs2gjAGBSIAwgsJwxBIodVVD5uQYYUAjA5EEYQCCZ3FB+DIEiZh6U8tMQM+kQgMmFMIDAMdnB/BgC2eJWDMedKwKhaHUKAwCXEAYQKCbTJ/W9KckubsVIvTPEsBWuSl0A4CbCAALDpDdJA28Vv2K0SUrOlGUxySeAyYkwgEnPGCOlO6TB9cWvHG+T6raTZVmVLwwAPIIwgEnNGCMNviOlO4tfuW47WYkplS8KADyGMIBJy9hZqf8tKdtf/MrJWbJizZUvCgA8iDCAScmZdfDN4nsMKCQ1MAUxgGAhDGBScdoHdEmD64pf2YrkBxNKVL4wAPAwwgAmDWNyUv9aKdNT/MqhmBMEwrHKFwYAHkcYwKRgcinntkCxIwpKUrjOGUMgxNsBQDDx6QffM+luZ9ZBmeJXjjZIydmMIQAg0AgD8C1jbGngHWmoq7QNxFql+umMIQAg8AgD8CWTG1LP2uVa8vQz+svy1yRJO83bQR//6EI1NzVOvIG67aR4G0EAAEQYgA+ZoV79y+WX6rs33al0ekhtbc5JvaOjQ60tTfp/37xMp5/0qbFXtiJO10GmHwaAYdwohW8YY2QG39H6VS/p29+9RZ/73Hl6/fXX1dHRoY0bN+rNN9/UUf90tM688Cv69f/8bvQGIkmpaR5BAAC2QhiALxg7K/WtllId6h8YlCTNnj1bX/va1zRjxgxNnTpVd955p+644w7tt99+uvb620duINGen36Yi2EAsDXLGFNCE2ygdkym3xlWOD+aYDab1bz3H6o3316n+fO213Gf+oReXbFKP334V3rjjTd077336tvfulqdf18mWWFnxsFoAe0IACCg+JoEzxpvtsFIJKIXf/tTrfz7Gu291x4Kh8P6ylXfVywWUzKZ1MaNG9Xc1CCFE063QQYSAoBtIgzAk4ydkwbekjJ9Yz7fPqVV7VNaJUn/+dBjuvq6W/TNb35TknTXXXfp2E99Umqcy/gBAFAAbhPAc0ymzxlEyM5MuOxDv/y1jj31CzrllFN0yy236Mgjj9QLL/xBf/7zXzRjxowaVAsA/seVAXiGsXPS4DvSUHdBy//i0ad03OmX6Pjjj9ctt9yiE044QUuWLNFjjz1GEACAInANFZ5ghnqlnpUFB4FHf/20Fp92sY4++mjdddddOu200/Szn/1MRx99tDo6OtTb21vdggFgEuE2AVxl7Kw0sK6omQb7+vo1e4+P6qP/+DHdf//9Ouuss/TjH/94xDL77LOPnnvuuUqXCwCTEmEArjDGOAFgYJ1kckWt++LLr2jvgxfrxRdf1F577aV0Oj3i+bvvvluf/exnlcvlFApx8QsAJkIYQM0ZO+M0EBynp8BENmzs1Ow9Pqq5c+cqkUiMer6rq0vRaFQrV65k7gEAKAANCFEzxhinTcDAO5Lskrczdc6ueuSRR/TII4/ItkdvJx6P6zOf+QxBAAAKxJUB1ITJDUkDb0vZgdI3YkWk5CxZ0WTlCgMAcGUA1eWMItiZH0WwjNwZa5bqpssKhStWGwDAQRhA1ZhcWup/W8oNlr6RUFSqnyEr2lC5wgAAIxAGUHHGGCm10flTztWAeKtUN02WxdUAAKgmwgAqymQHnZ4CuVTpGwnF8lcDaBsAALVAGEBFGGNLqQ1SqqO8DcWnSHVTmWAIAGqIMICyON0Fe6TU+oImFhpXOC7Vz5QVqatccQCAghAGUDKT6XcmFirnloAkJaZKiXbGBQAAlxAGUDSTSzldBUscQXBYOCElZ8oKjx5FEABQO4QBFMzYGWlwQ8EzC47PkuqmSfE2rgYAgAcQBjAhp3HgxnzjwDIHrIzUO20DwrGK1AYAKB9hAOManktgcH3RMwuOFpLqp0mxVq4GAIDHEAYwijO9cJ/TONAeKn+DkQYpOUNWKFr+tgAAFUcYwAgmO+iEgHImFNrMCkt120mxZq4GAICHEQYgKT+r4OB6KdNTga1ZzlDCialMLAQAPkAYCDhj55yRA9NdKrtxoCRFm5z5BGggCAC+QRgIKGNsZ2rh1EbJ2OVvMFIv1W3HCIIA4EOEgYBxQkC3lO4ob/jgzUIxp11AtIF2AQDgU4SBgDB21rkSkO6qQDdB5RsHTpNiLYQAAPA5wsAkZ3JpKdWZHzWwAm0CZEmJKfm5BJhZEAAmA8LAJGWyA86IgZneym001uJML8x4AQAwqRAGJpHhwYLSG6XsYOU2HG3I9xBgQiEAmIwIA5OAMbY0tMm5ElCJEQM3CyecHgLRZOW2CQDwHMKAjxk75zQITHdKJlu5DYeiTuPAaBONAwEgAAgDPmTsjHMVIN0tqQJjBGxmhaREe35qYRoHAkBQTKowYOyclEs5/eeNkWScE5wVci55h6K+/qZrsilnfIChTRXecig/fHA7wwcDQJ4xxjmf5FL5wdk298gKSaGIFE5Mms9MX4cBk8s4XeYyA1JusIBBdEIykYQUqZNiTVIk6flwYIyRsv3OlYBsf2U3bkWkRJszrfAk+YUGgFINf94O9TrnlGxKE119NaGoFK5zzivxFt/2trKMMZXofF4zTov5XqfvfLnd5kJRp898vFVWqPq5yKQ6nH0VcAne2DnnCsBQl5RLV7aQUMx53bFmbgcACDxnULZ8+6tyR2aNNkrxNt+NyuqrMGDSm6SBtZUZRndr8TapfnrVviGbVKc0uE6qnykr3jL2MsY4UwcPdUtDParMIEFbiNRL8Sm++yUFgGowds75XE53VX7jVkRKzpQVa6r8tqvAF2HA2Bmp7+0KTa+7DVZEaphV8YNn0t3SwNvOPyJJWY07jHzezjoBIN1d2a6Bm0UbpcQUWZH6ym8bAHzIDPVI/W9XtifWWKJNUnKG528feD4MmHSXEwQq2Wp+IrFmKTmrIlcJzFCv1L9m5IPNOznBI9vvJNJKjhI4zJLiLVJ8CtMJA0CesXNS/1vV/3I5QsgJBPHWGu6zOJ4NA8YYaXC988cNobjUtKOscOlpzmT6pL41GnW5P1Lv3Oqoxu0OK+z0DIi31aQdBAD4hbEzUs8qya5wO6xCJaY6o7l68DatJ8OAMcZpG5DqcLcQKyI1v6ekb9YmOyD1rlbF7/uPJxR12gPEW2gUCABbMbkhqef16t8WmEh8itM+zWOBwJthoH+tlNrodhmOUNQJBEXc7zG5lNS7Kt8vtcrCCWegoGij5365AMALnCsCr1fnamwpEu2y6qe7XcUIngsDJr1J6nvD7TJGiiSdWwYFnGxNLu1cEah6o5QGJ2FG6gkBADAOY4zU+3enp5aXNMyRFWt2u4phnrqpbOyM1P+m22WMtnnQn7r2bS5W/SBgOYMlJaYwgyAAFCLd4b0gIEn9b8lE6j3Ty8AzYcAYI/W9WZtL66UYWCsTbZAVGfskXNUgEI5LsVZnkCBGCgSAgphsShpY53YZYzO2EwgadvDE1V3PhAGlu6VMn9tVbFvfGpnm+aMOXHWCgOV0cYy3OuNfe+CXBQD8whjjzSvNW8r0OWPMeKDLoSfCwHA3Qq/LpZxbBtGG4YcqHgTCCecXI9Yky+IqAACUJNvvfGZ73eAGmViL61/4PBEGlOmvzsh71TC4cTgMVCwIWCHnKkCsddzbEACAInilR9pE7KFRXzLd4I0w4JeDJkmZXqe/qky++2CuvO1ZYal5J8YGAIAKMbkh79923lJqo+thwPUzkHPQqjEc70i9vf360MdO0F4HLdaeBx6tH93909I3NrihMkFAcrbhh0tZAOAX6c6a7eroz1yktnkLtfi0S0rfSKYv/yXTPa6PM+DMPVD9Rh65XE7p9JDq6+s0MDCo9x54tJ779X2a0tZS/MZCCal5rpQbci7xbP5788/FjjoYa5WVnFF8HQCAUcym12r2Jet/nnlOff0Duvv+h/XAnf9e+oaSs8ed0bYW3L9NkB2syW7C4bDq6+skSan0kHI5WyXnIDstWRFZ0aik5IinjDFOG4LhgJDZ4k82PwLWVvvNbJIx23GrAADGYfredLpZJ9q32djOGLumV1sPPvBDWrL0+fI3lB10Jpdziftnn0x5g0HYtq3d9jtSX/rX7454/ImnnlV8xgf0wC9+NfxY96Yevf8fjtWcPQ/RP194utqnlNqdw4z7y2ZZlqxQVFY0KSveKqtumqzkLFmNc2U1z5dadpWad5Yad5SSc6S66c4YAl4ZJhMAvCg7KHX9TdrwoszghvG/zOXKn4TItm3ttu8ifenKkd/0n3jqWcWn76UHfvFE2fsYxeWBkVy9MmDM+CfVQoVCIV3+hTP1+cu/rcsvOlOtLU16+S9/03FnXKqrrvi8Fh/18eFlW5qb9NJvf6p31nfomNMu1rGLDtF206aUtuNUh0w4Xlbtw6yQNNRTqymNAMB/svlzxeZQEHlDpnH70VcKKnC12TmvnKXPX/4tXf6FM9Xa0uycV07/oq76ykVafNShZe9jlFxKxhjXuhi62mbA2Bmp669lbyebzWqXfRfplOOP1Jknf0ofPvRkHXnYwbrxO1eMu87nLv2mPnrgviPCQlGMLQ1Vv+EjAGAbInXSFqHADLwjpTaUvdlsNqtdPnSETvn0UTrzpE/pw4ee5JxXrv3KqGWXLH1eN956X3ltBiSpZVfXpp539zZBhXJIJBLRZReeoRt+dK8O//T5+sCeu+n6b10+Ypl31neop9fpatLT26dnlr2oXebPLX2nfI0HAPdtdftAdgV6eSl/Xvn8Gbrhlp/o8E+fpw+8b3dd/+0vV2Tb43JxOH53rwzk0lL3qxXZVl/fgKbt+g+av+McLXv8x0om60c8/8JLr+isL3xdxhgZY3Tuacfpc2ccX/oO7Zy/+rECwGQXbZTizVK2Mp/NfX0DmrbLR5zzyhM/GXVekaRPHHuOXvzTcvUPDKqtpUkP3f197fOBBaXtsHknWZW6/Vwkl3sTVO7CxIWXXy1J2tjRrXB49DC+H3z/7vrjkgcqtj+JuQIAwBOijc6tgniLNPhOxcLAhV+6SpK0sXPs84okPf7TH1ZkX5Kc9mMucTcMhCrzwr/6rf/QI79+Rsse/7EOOeazuu0nP9P5Z366ItseVzgmqbG6+wAAOHKp0b2utggBmxvemQqdUL969Q3OeeWJn+iQo8/WbT9+SOefdUJFtj2uoIYBywrLhKJldau79Z4H9e833a0nf/YjvW/BLrronJN07Q136LOnHKNotIrzRDfMlhUjDABALZjuFdLAO84/xggBw8J1Ze/r3fPKrfnzysnOeeXUY6t3XglFXZ2czv1xBiKj78EU6rHfPKMLvnS17rn5au239/skSReefaJ6evt0z3/9d6UqHFuk/F84AIDDGCPT83eZ3tUy/W/LDK6XSXfJZHplsilJxgkBbXtI7XvKSrSO3Q2vzM/mx37zjC647Crdc/O3tN8+NTyvVCDElMMDYaC0/4AXXnpFx515qa75+sU6+oiPDT/e1NigC846Udd8/3blcpVpVTpKKOJa9w8AmIwsy3K+HOZSznw16S5navv+t6W+1XLCQJ001C0NbA4L3TKZARk7MzwIkRWKSFZpn88vvPR/Ou6ML+qaKy/R0Yu2OK80NeiCs0/UNd+/rXrnlTK+GFeC+3MTZPqlntfdLKF4sWZZjdu7XQUATComOyj1vVHi2pYUijrtuXKZijUirJnGHWVFkxMvVyUeuDJQ7xxAP3Fx/GgAmLTCiTLOB8aZDyYclxJtFS2r6qyI61cGXA8DlmVJiRKHBHZDKOrctwIAVIQxtszQJqlvTXnztCSmyEq0O5/Rlo++ZE4w+VIteOPGd7wt30rUB8P6Jaa4ftAAwO+cuWnSThuAoV5JZY6+l2iXlf9iaVmWTGKKNLiu7Dqrz5LipU6aVzmeCANWKCwTb5XSnW6XMgHLCS4AgJIYOydleqT0Jmc6+ErYIggMi7c6AxB5/UtmvFVWyL0uhZt5IgxIkuqnOQnRxbGZJ1Q/3RMHDQD8xBjjTNE7tCk/jHsFT9BjBQHlv2TWT5cG1lZuX5VmhaS6aW5XIclDYcAKRWWSs8toSVplkaS/2jYAgMuMncvfBthUXluA8SSmytpWY8F4mzTUI2X7K7/vSkjO9kw3ddcbEG7JijdLsRa3yxhDyBlxkLYCADAhkxtyphLuWSmlNlYpCLRvOwgo30A9OUseO9U5Yi2yYk1uVzHMG5FkS8mZToqrxi9PqRpmyQrH3K4CADzNZAedtl/VntE13jbmrYGxWOGYTHKW1L+mujUVIxSV6me4XcUIrg86NBaTG5I2rZRM1u1SpOTMgn/pACBojDHOyT/d6YweWDGbr8RudYqKtUh104q+UmtSHd5oP2BFpKZ5nvuC6ckwIOUDQc/r7l4hIAgAwJiMsZ22AOmuyn5Oh+NSrFmKNknpDmf7m0WbnIbcJd6ydT0QhKLOSIMeCwKSh8OAJBk7K/WsknKDNd6zJTXMcdowAACGGTvrnKAr2fvLCjkn+lizrEji3X3l0lLvKucf0QapfmbZbbecwY3eVM27HIbrpMYdPNNgcGueDgNS/hJUamPtBiWKJJ3Ggh5MbgDgFpNL50NAjyr2WRypz18FaJBljd3Iz/Suci6tJ2dVrBG3yQ1J/W/VqJeB5XQf9MAog9vi+TCwmcmlnTSXHajODqyQVD9z7PmxASCAhscHSHdV7sRpRZwAEGsq6EuXyfRLkbpxw0KpjDHO1Y3+tSp79MPxROqdEBOOV2f7FeSbMCBtbqjSI6Uq2FrVijiTWiSmePbyDQDUkvNZ25tvFFihUQLDdVKiVYo0eOoLl3Pbo9M5r1Sq0frmcWmijZ56rdviqzCwJZMbyh/ArtIOYLTRCQE+OlgAUE3G5JxhgtMlfq6OJdroDLkbqavM9qrk3V4RHaV92bQizhDI8VZf3mb2bRjYkrEzUnZQyqacv+2MJNu5rWVZkhWWIgknmUbqpHCcAAAAecPfjtObVJlL5pbTBTDe4ssT4/AkSrn8eSWXkkzOaTBpWZJCUiiSP5/USZGErJKnXvaGSREGAADFM8Z+9xJ5JRoFWmHn23GshXlcfIab5AAQME7juU1OTy2TK3+DoVj+tmsTV119ijAAAAExfF88tVGyh8rfYKTemQwoUk8I8DnCAAAEgMkOSoPrKzNkcLRJSrTKCicmXha+QBgAgEnM5NLOlYCyu2OHpHiL0yjQ543lMBphAAAmIWNnnRAwtKm8DYWi+UaBzRUf+AfeQRgAgEnEGNvpHZAus4dAOOGEAMZiCQTCAABMAsPD66Y6yushEIpLdVNpFBgwhAEA8LF3ewhsKG8qYSsi1bXTPTCgCAMA4FMmOyANbiivh4AVkuJT8pO00SYgqAgDAOAzJpd2QkBZMwla+d4BUxgtEIQBAPALY3LS4EanbUA5ok1SXTtdBDGMMAAAPmAyfdLAO+XNJhipl+qmMlgQRiEMAICHGTvrjByY6S19I+G4lJgqK5qsXGGYVAgDAOBBTi+BXicIlNpVMBSVEu2MFYAJEQYAwGOMnXFuCZTaQJAeAigSYQAAPGJ4auHBDZLsErZgOaMGxtvoIYCiEAYAwANMbkgaWCflBkvbQKxJStBDAKUhDACAi4wxzjwCqQ6VNJdAuE6q305WOF7x2hAchAEAcInJpqTBdVIuXcLaIWcOgVgzjQNRNsIAANSYM7NgR35mwRJEks7VAG4JoEIIAwBQQyY74LQNKGVSISss1U2jqyAqjjAAADVQ9lDC0UapbpqsEB/bqDx+qwCgBJlMRq+88or+9Kc/KZ1O6/3vf7/23nvvMZctayhhK+LcEog2lFkxMD7CAAAU6eKLL9bNN9+sdHpkw79nn31WCxcuHP63MbYTAjI9pe0o1uJMKGQxZgCqi6GpAKAIa9eu1fe+9z2dd955Wrp0qfr7+2WMUX19vf7whz8ML2dyaal3dWlBIBSVGubIqt+OIICaIAwAQBFSqZQkqbGxUY899piuv/76UcuYoU1OELCHit9BvE1qnCsrUl9uqUDBuE0AAEXYYYcdtM8+++gb3/iGLMvSkUceOfycMUZmYJ0zpHCxwnGpbrqsCNMLo/a4MgAARQiFQlq2bJm6u7v1yU9+cuST6c4SgoDlzCzYsANBAK7hygAAFCkcDqu5uXn0E8X2FgjXSfXTZYVjlSkMKBFhAABqjqGE4S2EAQAogbEzUnZQijQWt2K4TkrOYChheAptBgCgSJ3r1+iOm6/TqtVrRjz+zO/+oJ8/8htlMuMMNRxvdboMEgTgMYQBACiQMUZmcIOOPvoYnfOFr+q1v7+hWMy5359IJPTfT/xWx57yeV17/W0jV7RCUnKmrLpp3BaAJxEGAKAAxs5K/WukdKf++Kfl+va3v61UKqX/+q//kiR1dHQolUrpkEMO0Ut/Xv7uiuG401MgWuTtBKCGLGOMcbsIAPAykxmQBt6WTE6StN/HjtOK19eosXH0CX79+vX64gWn65tXXJQfTniqLIvvXfA2wgAAjMMY44wdkNo44vHVb7yle3/63xoYTI1aZ9aM7XTycUeqYdp7ZMWaalUqUBbCAACMwdg5aWCtlO0vbsVQzGkfEI5XpzCgCuhaCABbMdlBqf/t4gcRijY50w1zWwA+QxgAgC2YoR7nikBRLKluGoMIwbcIAwCgze0DuqTUhuJWDEWl+pnMKwBfIwwACDxjjDS4XhrqLm7FaIMzt4AVrkpdQK0QBgAEmjG2c1sg01fciompUryV2wKYFAgDAALL2Dmp/y0pN1j4SlbE6S0QqateYUCNEQYABJKxM1Lfm5I9VPhKkaRzWyDERycmF36jAQSOyaWkvreK6zqYmCLFp3BbAJMSYQBAoJhMvzOGgOzCV6qfLivWXLWaALcRBgAERvFjCFhScpasaLJqNQFeQBgAMOmVNIaAFZaSsxk/AIFAGAAwqZU0hkAo6gSBcKxqdQFeQhgAMGmVNIZAOOHcGqDHAAKE33YAk1JJYwhEGqTkDCYaQuAQBgBMOiWNIRBrkeqm0XUQgUQYADCplDaGQLsUbyMIILAIAwAmDcYQAEpDGAAwKZhMn9NGoGCMIQBsRhgA4HvvXhEoEGMIACMQBgD4mskM5K8ImMJWYAwBYBTCAADfMtkBqf9NFRwEGEMAGBPvCAC+ZLKDTq+BQoMAYwgA4yIMAPAdk03lrwgU2GuAMQSAbSIMAPAVk0s7QcAUGATibVKinSAAbAPXywD4hskNSX1rJJMrbIV4K0EAKABhAIAvFB0EYi1SYipBACgAYQCA5zlzDawpfIjhWDNtBIAiEAYAeFrRQSDaJNVtRxAAikAYAOBZxs7mZx/MFLZCtNGZa4AgABSFMADAk5wgsKbwaYijDVL9DIIAUALCAADPMXbO6T5YaBCIJAkCQBkIAwA8xZh8EMilC1shUi8lZzKyIFAG3j0APMMY2xliOJcqbIVwnTPXAEEAKAvvIACeYIztzD6YGyxshXBCaphNEAAqgHcRANc5QeBtKTtQ2ArhuDMNMUEAqAjeSQBcZYyRBt+Rsv2FrRCKSck5skLh6hYGBAhhAIC70l3SUE9hy4ZiUgNBAKg0wgAA15hMn5TaUNjCoWg+CDDZKlBphAEArjC5lNNOoBAEAaCqCAMAas4ZXfAtSWbiha2I01gwFK16XUBQEQYA1NRwF8JCJh6yws4VgXCs+oUBAUYYAFAzxhhp4J3CBxVKziQIADVAGABQO+lOKVNgz4G66bIi9dWtB4AkwgCAGjFDvVJqY2ELx1tlxZurWxCAYYQBAFVnsilpYG1hC0eSUmJqdQsCMAJhAEBVGTvrNBgspOdAKJafgZCpiIFaIgwAqJqiew4w3wDgCt51AKrC6TmwrsCeA1a+5wBjCQBuIAwAqI50h5TpLWzZ+u3oOQC4iDAAoOKcngMdhS0cb5MVo+cA4CbCAICKKq7nQIOUaK9uQQAmRBgAUDHGzhTRcyAuJWfQcwDwAMIAgIoofs6BWfQcADyCdyKAsjk9B9ZKuXQBS1tSchazEAIeQhgAUL50p5TpK2zZ+umyInXVrQdAUQgDAMpisgNFzDnQJivWVN2CABSNMACgZM5QwwX2HIjScwDwKsIAgJIMjzBYSIPBcFyqp+cA4FWEAQClSXdK2f6Jl7PCToNBeg4AnsW7E0DRimonQM8BwPMIAwCKYuxc4e0EElPpOQD4AGEAQMGGxxMopJ1AJCnFW6tfFICyRdwuAEDtGWNL2QHJzkrGdh60QlIoIkXqx7+/n+4qsJ1AhAaDgI8QBoAAMLm0M4vgUJ8zrXB2cNvLR+qkaKMUa5ASU2SF4zLZQSm1obAdJmfKCoUrUDmAWrCMMQXMKALAb4wxzjf5/nVOy39JkqWCJhHaetl4m5RozY8yaG97tcRUWYm2kmoG4A7CADDJGGOkwfVSz2rJHqrsxkNxqa5VMrmxn48k890IuT0A+AlhAJhETDYlda+QhjZVd0eJVqd9wZasiNQ4l9sDgA8RBoBJwBjjdPfrWaUJL+NXSigmJVo0fCuhYXu6EQI+RRgAfM4Y41wNGFzvwt4tqa5dqmuXlZjiwv4BVAK9CQAfM8aWuv4qpTonXrg6FUiDGySTk4m3MuQw4FO8cwGfMsZIXa+6GAS2kOqUul4VFxoBfyIMAH7V92bh8wPUQmqjUxMA3yEMAD5khvqk3tVulzFa72qZTJ/bVQAoEmEA8BmnncDf3C5jfJ1/c2oE4BuEAcBveldLuW0PJ+yq3KDU+4bbVQAoAmEA8BFjZ6W+t90uY2J9bzu1AvAFwgDgJwPvqPC5Bdxk52sF4AeEAcAnjDH+uCqwWd/bdDUEfIIwAPhFukuy0zXb3cBgSnM/dKwu/caNpW3ATjs1A/A8wgDgF+luOdMK18ZV379b++61exlbsPI1A/A6wgDgF0O9qlV7gRWvr9HfVr6hwz66XxlbMfmaAXgdYQDwAWOMVOZgPrZta7ePnKQvXXXziMefWPJ7xecerAd++T/Dj/3zN2/S1Zd/tqz9SZIy/bQbAHyAMAD4QXZA5V4VCIVCuvyCk/WDu3+urm7nG/vL//eajjvna7rq8s9q8aKDJUm/eOIZ7TRvtnZ+z/blVi3JztcOwMuYwhjwAZPqkDqXl72dbDarXQ48UacsPkxnfvpwfXjRuTry0AN049WXDC/z5W/9QD958NcKh0Pq6x9UJpvVJeccr69dfHppO23bjemNAY8jDAA+YAY3VGwI4h/e8wtd8e1bNGtGu+bOnqGHbrtK4XB4zGXv/M9H9Ze//V3/72vnl77D1l1k1U0tfX0AVcdtAsAPKpjZTzr6EA2kUjJGuvemr48bBCqG7xuA50XcLgBAAazK5fYLr/iuJGljZ7fCoW1v97TjP1n+DitYO4Dq4F0K+IFVmW/vX/3OrXrkyWVa9vAPlc3mdNv9j1Rku9tUodoBVA9hAPCDaLLsTdx67y/17z+8Xw/feY3et8d8XXTWcbr2pnuVyVR5QqEK1A6guggDgA9Y4ZgUipa8/mNP/a8uuOK7uueGr2q/D+4hSbrwjGPU09evex58olJljhaKOrUD8DTCAOAX0caSVnvhT3/Tced8Tddc8Tkd/cl/GH68qTGpC04/Rtf8x4+Vy+UqVeVIsabqbBdARdG1EPAJ07tG6l3tdhnFadxBVuMct6sAMAGuDAB+UdfudgXF82PNQAARBgCfsCJ1UrzV7TIKF291agbgeYQBwE+SM92uoHB+qhUIOMIA4CfxFimccLuKiYUTTq0AfIEwAPiIZVlS8zy3y5hY83ucWgH4AmEA8Bkr0SbVb+d2GeOrny4r4aO2DQAIA4AvNc2TwnG3qxgtHJeadnS7CgBFIgwAPmSFwlLrLpK8dCnecqYrDjEXAeA3hAHAp6xYk9S2q9tlvKttV6cmAL5DGAB8zEpMkVp3lbtXCCypdTenFgC+xHDEwCRg0pukzlckU6U5BsZjhaW23WXFm2u7XwAVRRgAJgmTS0vdK6V0Z212GG+TWt4jy4sNGQEUhTAATCLGGCnVIXW/JplsdXZiRaSW+bKYdwCYNAgDwCRk7IzUv9b5Y2cqs9FQVErOkJIzZYUildkmAE8gDACTmDFGSnc5oSDdtcUzlqTx3vpbPRdvdUJAvJVRBYFJijAABISxc1K2X8r0S5k+aagvf9Vg80eA5Xz7jzVI0QYpmpQiScYNAAKAMAAAQMAxzgAAAAFHGAAAIOAIAwAABBxhAACAgCMMAAAQcIQBAAACjjAAAEDAEQYAAAg4wgAAAAFHGAAAIOAIAwAABBxhAACAgCMMAAAQcIQBAAACjjAAAEDARdwuYIRUh2RsScb5e/hnM8Zjtszmx2Ukk9tqOfPuuia3xc+2JPvdx7fetrFHrGPMVs/JSPYW25CR7NxW2zDj/LHz5W5R46h/T/DH3mrZXO7dx3K2U9sYz5ktn8vlX/eW/7bNVs/lX7tt8i/byNibH8//OzfyZ+e/OP93ztmeMUbavJyd34ZR/m8jO2fL2LZkG9m53PDjZst/23Z+ufxryT+3eXvDPxsjYzuv1djOz3bOzj/u7MvO/2yMkZ179+fc5nW2WNYYI3uLn99dz5a9+WfbKP/bICMpm//blpQb5/HNz2293OZlc1std6Ux1Xi31R7vb97fvL89+/7mygAAAAFHGAAAIOAIAwAABBxhAACAgCMMAAAQcIQBAAACjjAAAEDAEQYAAAg4wgAAAAFHGAAAIOAIAwAABBxhAACAgCMMAAAQcIQBAAACjjAAAEDAEQYAAAg4wgAAAAFHGAAAIOAIAwAABBxhAACAgCMMAAAQcIQBAAACjjAAAEDAEQYAAAg4wgAAAEFnPCKVSpmvf/3rJpVKuV1K4HEsvGOyHIvJ8Domw2swhtfhJV56DZYxxrgdSCSpp6dHzc3N2rRpk5qamtwuJ9A4Ft4xWY7FZHgdk+E1SLwOL/HSa+A2AQAAAUcYAAAg4AgDAAAEnGfCQDwe19e//nXF43G3Swk8joV3TJZjMRlex2R4DRKvw0u89Bo804AQAAC4wzNXBgAAgDsIAwAABBxhAACAgHMtDFx66aU68MADddJJJ2loaGj48Ww2q9NOO00HHnigLrroIrfKC5TxjsXKlSu11157KZFIqK+vz8UKg2O8Y/Hoo49q4cKFOuCAA3TBBRe4WOH4in1Pf/e739X++++vI444Qps2bXKj5DEV+zoaGxt10EEH6aCDDtKf//xnN0oepdj3tN+OxXivw4vHQir+fe3G8XAlDPzxj3/UunXr9Mwzz2j33XfXT3/60+HnfvnLX2r27Nl65plnNDAwoN/97ndulBgY2zoWM2bM0JIlS7Tffvu5WGFwbOtYLFiwQE8//bSWLl2qzs5OPf/88y5WOlqx7+kNGzbol7/8pZYuXaoTTjhBN954o4vVv6uUz6ZddtlFS5Ys0ZIlS/Te977XrdKHFfue9uOxGO+zyWvHQir+fe3W8XAlDCxbtkwf//jHJUmf+MQnRpzwt/UcKm9b/9/19fVqbm52q7TA2dax2H777RWJRCRJ0Wh0+GevKPY9/fzzz+uggw6SZVmeep+X8tm0cuVKfeQjH9HnPvc5pVKp2he9lWLf0348FuN9NnntWEjFv6/dOh6uhIHu7u7hcZibm5vV2dlZ0HOoPP6/vaOQY/HCCy9o48aN2muvvWpd3jYV+5726u9dKZ9Nr732mp5++mnNmDFDN910U+2L3kqx/7d+PBbj8dqxkIp/X7t1PFwJA62trerp6ZHk/Ee1tbUV9Bwqj/9v75joWLz55pu66KKLdOedd7pQ3bYV+5726u9dKZ9NU6ZMkSQtXrxYL730Um0LHkOx/7d+PBbj8dqxkIp/X7t1PFwJA/vtt59+9atfSZKeeOIJ7b///gU9h8rj/9s7tnUs+vr6dOKJJ+oHP/iBpk6d6laJ4yr2Pb333ntryZIlYy7vpmJfR39/v3K5nCTp6aef1vz582tf9FaKfU/78ViMxYvHQir+fe3W8XAlDOy1116aPn26DjzwQL3yyis65phjdM4550iSFi1apDVr1ujAAw9UXV2dPvzhD7tRYmBs61h0dXXpYx/7mF5++WUtWrRIjz32mMvVTm7bOhY33HCDVq5cqQsuuEAHHXSQfvvb37pc7UjFvqenTp2qRYsWaf/999d9992n8847z+VX4Cj2daxYsUL77LOPPvKRj+jRRx/1RA+oYt/TfjwWY70OLx4Lqfj3tVvHg+GIAQAIOAYdAgAg4AgDAAAEHGEAAICAIwwAABBwhIEAuPPOO9XS0lKRba1atUqWZSkSieitt94a8dzatWsViURkWZZWrVo14rkHH3xQBx10kJqbm9XQ0KA999xT3/jGN4YH1KhkjUDQnHbaabIsS+eee+6o58477zxZlqXTTjtt+LF169bpwgsv1Lx58xSPxzVnzhwtWrRITz755PAyc+fO1fe+970aVA8vIAygJDNnztTdd9894rG77rpLs2bNGrXsFVdcoeOPP1777LOPHnvsMf3lL3/Rddddp5dffln33HNPrUoGJrU5c+bo/vvv1+Dg4PBjqVRK9913n7bffvvhx1atWqUPfvCDeuqpp/Sd73xHf/7zn/X444/r4IMP1vnnn+9G6fAAwoAPPP744zrggAPU0tKiKVOm6IgjjtDKlSslSUuWLJFlWeru7h5e/qWXXhr+dr5kyRKdfvrp2rRpkyzLkmVZuvLKKyU5fXVPOeUUtba2qr6+XocddphWrFhRUE2nnnqq7rjjjhGP3XnnnTr11FNHPPbcc8/p6quv1nXXXadrr71WCxcu1Ny5c3XIIYfowQcfHLU8gNJ84AMf0Pbbb6+HHnpo+LGHHnpIc+bMGTF89eYrBc8995yOPfZY7bzzztpjjz10ySWX6H//93/dKB0eQBjwgf7+fl1yySV6/vnn9eSTTyoUCulTn/qUbNuecN2FCxfqe9/7npqamrR27VqtXbtWl156qSTn0uIf/vAHPfzww1q2bJmMMfrkJz+pTCYz4XaPPPJIdXV1aenSpZI0POvWokWLRiz3k5/8RA0NDeMOnMGtAaByTj/99BEh/fbbb9cZZ5wx/O/Ozk49/vjjOv/885VMJketz/sxuLw19RnGdMwxx4z492233aZp06bplVdemXDdWCym5uZmWZal6dOnDz++YsUKPfzww3r22We1cOFCSc6Je86cOfr5z3+uxYsXb3O70WhUJ598sm6//XYdcMABuv3223XyyScrGo2OWG7FihWaN2/eqMcBVN5nPvMZffnLXx5u2/Pss8/q/vvvHx7e9rXXXpMxRrvuuqu7hcJzuDLgAytXrtSJJ56oefPmqampSTvuuKMk6Y033ih5m8uXL1ckEtG+++47/NiUKVO0yy67aPny5ZKkww47TA0NDWpoaNAee+wxahtnnnmmHnjgAa1bt04PPPDAiG8gmxljZFlWyXUCKFx7e7sOP/xw3XXXXbrjjjt0+OGHq729ffj5zQPO8p7E1rgy4AOLFi3SnDlz9KMf/UgzZ86UbdtasGCBhoaG1NDQIOndN7mkgi7zjzcK9ZYn71tvvXW4MdJY3+wXLFigXXfdVSeccIJ22203LViwYNRMYTvvvLOWLl2qTCbD1QGgBs444wxdcMEFkqQbb7xxxHM77bSTLMvS8uXL9U//9E8uVAev4sqAx3V0dGj58uX6yle+on/8x3/Ubrvtpq6uruHnN890tXbt2uHHtj4hx2Kx4dm8Ntt9992VzWb1+9//fsS+Xn31Ve22226SpFmzZmn+/PmaP3++dthhhzHrO+OMM7RkyZIxrwpI0oknnqi+vr5x5xbfsuEjgPJ94hOf0NDQkIaGhnTooYeOeK6trU2HHnqobrzxRvX3949al/djcBEGPK61tVVTpkzRLbfcotdee01PPfWULrnkkuHn58+frzlz5ujKK6/Uq6++qkceeUTXXXfdiG3MnTtXfX19evLJJ7Vx40YNDAxop5120lFHHaWzzz5bS5cu1csvv6yTTz5Zs2bN0lFHHVVwfWeffbY2bNigs846a8zn9913X1122WX64he/qMsuu0zLli3T6tWr9eSTT2rx4sW66667SvuPATCmcDis5cuXa/ny5QqHw6Oev+mmm5TL5fShD31IDz74oFasWKHly5fr+uuvZ5bYACMMeFwoFNL999+vF154QQsWLNDFF1+sa6+9dvj5aDSq++67T3/961/1vve9T9dcc43+7d/+bcQ2Fi5cqHPPPVfHH3+8pk6dqu985zuSpDvuuEMf/OAHdcQRR+jDH/6wjDF69NFHi7qcH4lE1N7erkhk/DtO11xzje699179/ve/16GHHjrcjWnPPfekayFQBU1NTWpqahrzuR133FEvvviiDj74YH3xi1/UggULdMghh+jJJ5/UzTffXONK4RVMYQwAQMBxZQAAgIAjDAAAEHCEAQAAAo4wAABAwBEGAAAIOMIAAAABRxgAACDgCAMAAAQcYQAAgIAjDAAAEHCEAQAAAu7/AwkOO1wbylQoAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "tp.plot_graph(\n",
    "    graph = boot_graph,\n",
    "    val_matrix= val_mat,\n",
    "    link_width = boot_linkfreq,\n",
    "    var_names=dataframe.var_names,\n",
    "    vmin_edges=0.,\n",
    "    vmax_edges = 0.2,\n",
    "    edge_ticks=0.05,\n",
    "    cmap_edges='OrRd',\n",
    "    vmin_nodes=0,\n",
    "    vmax_nodes=.2,\n",
    "    node_ticks=.1,\n",
    "    cmap_nodes='OrRd',\n",
    "    ); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
